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FOREWORD 


The  theoretical  basis  of  the  TIGER  code  for  calculating  the  thermo- 
dynamic state  in  a nonideal  heterogeneous  mixture  was  formulated  by 
Dr.  S.  R.  Brinkley,  Jr.  The  original  version  of  the  code  was  developed 
and  documented  from  July  1966  to  November  1968  at  Stanford  Research 
Institute  by  W.  H.  Zwisler,  W.  E.  Wiebenson,  and  L.  B.  Seely  for  Ballistic 
Research  Laboratories  under  Contract  No.  DA-04-200rAMC-3226(X)  monitored 
by  Dr.  S.  M.  Taylor. 

Development  of  the  code  after  1968  was  supported  by  Lawrence  Livermore 
Laboratory  under  Contract  AT(04-3) -115,  Agreement  No,  89,  P.O,  No.  5411209, 
with  Mr.  M.  Finger  and  Dr.  E.  Lee  as  Technical  Monitors,  New  routines 
were  also  added  to  the  code  as  a result  of  work  performed  for  Picatinny 
Arsenal  under  Contract  DAAA21-71-C-0454  monitored  by  Mr.  J.  Hershkowitz 
and  work  performed  for  Naval  Ordnance  Laboratory  under  Contract  N60921-72- 
C-0013  monitored  by  Dr.  S.  J.  Jacobs,  The  present  documentation  of  TIGER 
was  written  for  Lawrence  Livermore  Laboratory,  for  Naval  Ordnance  Labora- 
tory, and  for  Picatinny  Arsenal  under  the  contracts  cited  above  to 
provide  an  up  to  date  version  of  the  code. 
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GENERAL  ABSTRACT 


TIGER  is  ft  digital  computer  program  written  in  FORTRAN  TV  for  calcu- 
lating the  thermodynamic  state  attained  in  a heterogeneous  system  of  known 
atomic  composition  containing  gases,  liquids,  and  solids  with  arbitrary 
equations  of  state.  As  currently  arranged,  the  program  is  modular, 
includes  52  routines,  and  can  be  applied  to  systems  containing  up  to  30 
gaseous  and  10  condensed  constituents  composed  of  up  to  10  chemical 
elements.  The  memory  required  to  perform  a calculation  depends  on  the 
computer.  For  example,  on  the  CDC  6400  computer,  40K  words  are  required. 
The  input  instructions  are  written  in  a quasi-f ree-f ield  format  to 
simplify  the  input  of  data  into  the  computer. 

The  TIGER  documentation  consists  of  four  volumes.  Volume  I presents 
the  theoretical  basis  of  the  code  and  its  application  to  the  calculation 
of  the  detonation  parameters  of  condensed  explosives.  The  equations  used 
by  Brinkley  to  calculate  the  tnermodynamic  state  in  a nonideal  hetero- 
geneous system  in  chemical  equilibrium  are  derived  and  extended  to  treat 
systems  in  partial  equilibrium.  A brief  discussion  of  the  Chapman-Jouguet 
(CJ)  theory  of  detonation  is  followed  by  an  account  of  the  methods  used 
to  calculate  conditions  in  the  CJ  wave  and  the  properties  of  the  detona- 
tion products  along  Hugoniot  curves  and  isentropes.  Volume  II  presents 
a summary  of  the  formulas  and  relationships  used  in  TIGER.  Volume  HI 
presents  the  FORTRAN  code  and  flow  charts  needed  to  understand  the  pro- 
gram in  detail,  and  Volume  IV  is  a ..ser's  guide  that  explains  how  to 
prepare  input  cards  and  interpret  the  output  of  a calculation. 


VOLUME  r 


ABSTRACT 


Volume  I contains  the  theoretical  basis  of  the  TIGER  code  and  its 
application  to  the  calculation  of  the  properties  of  detonating  condensed 
explosives.  The  equations  used  by  Brinkley  to  calculate  the  thermodynamic 
state  in  a nonideal  heterogeneous  system  of  known  atomic  composition  in 
chemical  equilibrium  are  derived  and  extended  to  treat  such  systems  in 
partial  equilibrium.  A brief  discussion  of  the  Chapman-Jouguet  (CJ) 
theory  of  detonation  is  followed  by  an  account  of  the  methods  used  to 
compute  conditions  in  the  CJ  wave  and  the  properties  of  the  detonation 
products  along  Hugoniot  curves  and  isent ropes.  The  presentation  in  this 
volume  will  be  of  interest  primarily  to  the  reader  concerned  with  under- 
standing the  theoretical  background  of  the  code  rather  than  to  the  user 
of  the  program. 
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I -A . INTRODUCTION 
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TIGER  is  a digital  computer  code  written  in  FORTRAN  IV  for  calcu- 
lating the  thermodynamic  state  of  a nonideal  heterogeneous  system  of 
known  composition.  It  was  developed  at  Stanford  Research  Institute 
specifically  for  detonation  calculations  after  experience  with  the  RUBY 
code  at  Lawrence  Livermore  Laboratory  and  the  BKW  code  at  Los  Alamos 
Scientific  Laboratory  had  led  to  the  conclusion  that  a more  versatile 
code  was  required  to  perform  routine  and  research  calculations  on  con- 
densed explosives.  While  RUBY  is  limited  fay  its  inability  to  treat 
certain  explosive  compositions,  both  RUBY  and  BKW  are  restricted  by  the 
inflexibility  of  their  interlocking  subroutines  which,  for  example,  pre- 
vent use  of  a new  equation  of  state  in  a calculation  without  complete 
reprogramming. 

The  TIGER  code  .as  constructed  to  avoid  the  problems  associated  with 
interdependent  subroutines.  The  program  was  written  in  modular  form  so 
that  the  thermodynamics  used  to  calculate  the  state,  the  hydrodynamics 
used  to  calculate  detonation  parameters,  and  the  equations  of  state  used 
to  describe  the  properties  of  the  system  are  treated  separately  in  dif- 
ferent parts  of  the  code.  Because  of  this  separation,  TIGER  can  be  best 
described  as  a general  code  for  calculating  the  thermodynamic  properties 
of  a nonideal  heterogeneous  mixture,  described  by  an  arbitrary  equation 
of  state,  with  the  capability  of  calculating  detonation  parameters  pro- 
vided by  the  hydrodynamic  option.  Whereas  the  hydrodynamic  problems  are 
computationally  trivial  and  are  solved  in  the  executive  part  of  the  pro- 
gram, the  thermodynamic  problems  are  complex  and  are  solved  in  the  sub- 
routines THERU0  and  ECSBIPG,  which  constitute  the  largest  part  of  the 
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program , The  gaseous  and  condensed  equation 
named  STATE  G and  STATE  C.  They  are  called  on  by  THERM0  and  ECOMP0  in  a 
thermodynamic  calculation  whenever  equation  of  st?te  data  are  required. 
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I -B.  THEORETICAL  BASIS  OF  TIGER 
FOR  EQUILIBRIUM  CALCULATIONS 


1 . Introduction 

This  section  of  the  TIGER  documentation  presents  the  theoretical 
basis  of  the  methods  developed  by  R.  S.  Brinkley,  Jr.  to  calculate  the 
thermodynamic  properties  of  a nonideal  heterogeneous  system  of  known  atomic 
composition  containing  gaseous,  liquid,  and  solid  phases.  The  thermodynamics 
is  formulated  with  an  arbitrary  equation  of  state  to  make  TIGER  applicable 
to  a wide  variety  of  problems,  and  the  most  suitable  equation  of  state 
available  should  be  used  to  perform  the  calculations  in  a particular 
application. 

The  heterogeneous  system  is  assumed  to  be  in  mechanical  and  thermal 
equilibrium,  but  not  necessarily  in  chemical  equilibrium.  Since  pressure 
and  temperature  are  well  defined  in  this  case,  the  thermodynamic  state 
of  the  mixture  can  be  defined  by  the  values  of  the  state  variables, 
pressure  p,  temperature  T,  volume  V,  entropy  S,  internal  energy  E,  and 
by  the  mole  numbers  of  the  constituents.  Other  thermodynamic  variables 
such  as  the  enthalpy  H and  the  Gibbs  free  energy  G,  can  then  be  calculated 
with  thermodynamic  identities.  The  problem  addressed  by  Tiger  is  that  of 
computing  the  thermodynamic  state  when  the  gross  composition  and  a com- 
plete equation  of  state  of  the  mixture  are  known.  An  assumption  about 
the  attainment  of  chemical  equilibrium  and  the  specification  of  two 
independent  state  variables  are  required  to  solve  this  problem.  Section 
I-B  presents  methods  used  to  calculate  the  thermodynamic  jtate  when  the 
composition  is  assumed  to  be  in  a state  of  chemical  equilibrium  and  the 
pair  of  variables  can  be  chosen  from  the  following  set:  [(p,T) , (p,h) , 

(p,s) , (v,T)  , (v,e) , (p,v) , and  (v,s)],  where  h,  s,  v,  and  e denote  the 
specific  values  of  the  enthalpy,  entropy,  volume,  and  internal  energy. 
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Section  I-C  presents  the  extension  of  these  methods  to  nonequil ibrium 
systems  that  are  partially  frozen  and  partially  in  chemical  equilibrium, 
and  Section  I-D  presents  the  method  used  to  calculate  the  thermodynamic 
state  in  systems  that  are  completely  frozen  in  metastable  equilibrium. 

Assumptions  about  the  composition  are  required  to  calculate  the 
thermodynamic  state  when  kinetic  processes  are  not  treated  explicitly. 

The  assumptions  used  in  TIGER  are  related  implicitly  to  rate  processes 
and  cover  a wide  spectrum  of  chemical  kinetics.  As  a consequence,  the 
results  based  on  them  provide  a means  of  modeling  chemical  kinetic 
processes  with  the  TIGER  code.  When  the  system  is  assumed  to  be  in 
chemical  equilibrium,  the  composition  is  unknown  and  must  be  determined 
in  the  calculations  of  the  thermodynamic  state.  The  equilibrium  compo- 
sition is  calculated  with  the  equilibrium  conditions  and  the  stoichio- 
metric conditions  that  express  the  conservation  of  mass  for  the  system 
in  terms  of  its  molecular  and  atomic  composition.  When  the  system  is 
assumed  to  be  frozen,  the  composition  must  be  chosen  to  satisfy  the 
Stoichiometric  conditions;  when  it  is  assumed  to  be  in  partial  equili- 
brium, the  frozen  composition  is  chosen  and  the  remainder  is  calculated 
with  the  equilibrium  and  stoichiometric  conditions  as  in  the  previous 
case. 

2 . Description  of  the  Composition  ir.  Terms  of  Components 

It  is  convenient  for  computational  purposes  to  formulate  a general 
method  for  describing  the  composition  of  a heterogeneous  system  in  terms 
of  its  gross  composition.  The  possible  species  that  make  up  the  system 
are  restricted  to  a set  chosen  on  the  basis  of  chcii  ical  intuition  and  from 
the  results  of  previous  calculations.  We  will  consider  the  set  of  t 
possible  species,  s gaseous  and  (t  - s)  condensed,  Formed  from  c 
different  chemical  elements.  These  species  will  be  distinguished  wi th  a 
constituent  index  i and  the  elements  with  an  atomic  irdex  k.  Thus, 
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the  gaseous  constituents  are  labelled  with  i - 1,  . ..,  s,  the  condensed 
constituents  with  i = s + 1,  and  the  elements  with  k = 1,  ....  c. 

Condensed  species  are  assumed  to  be  present  as  pure  phases  to  exclude  the 
consideration  of  solid  and  liquid  solutions.  The  initial  composition  in 
a mass  of  mixture  Mo  is  described  by  parameters  specifying  the  gross 
composition  of  the  system.  The  final  composition  obtained  by  computing 
the  thermodynamic  state  of  the  system  is  defined  by  the  mole  numbers 
n (i  =1,  . . . , s)  of  the  gaseous  constituents  and  by  the  mole  numbers 
n*  (i  - s + 1 , . . . , t)  of  the  condensed  species.  The  phase  rule  imposes 
a restriction  on  the  number  of  the  variables  n*  that  may  be  nonzero. 

X 

When  all  the  n,  are  zero,  the  system  is  homogeneous  and  consists  of  a 
gaseous  phase  only. 

The  parameters  used  to  describe  the  gross  composition  of  the  mixture 
are  related  to  the  mole  numbers  of  the  constituents  by  the  law  of  con- 
servation of  mass,  and  such  a relationship  is  used  in  the  computation  of 
the  mole  numbers  and  n^  to  ensure  that  mass  is  conserved  in  our 

closed  system.  Since  iterative  procedures  are  in  general  necessary  to 
compute  the  composition,  it  is  important  to  express  the  conservation  of 
mass  relationship  in  terms  of  the  most  suitable  representation  of  the 
gross  composition  for  performing  the  iterations.  It  has  been  found  from 
experience  that  calculation  of  and  n*  usually  proceeds  most  readily 

when  the  gross  composition  is  expressed  in  terms  of  the  most  abundant 
species  in  the  system. 

Let  M and  Q denote  the  mass  and  number  of  gram  atoms  of  element 
k k 

k present  in  a mass  Mo  of  our  mixture.  Then  Mo  and  the  c values 
of  are  related  by  the  equation 

Mo  ■ ^ \ lt-0-1) 
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expressing  the  conservation  of  mass.  The  gross  composition  of  the  system 


can  be  described  by  using  the  c values  of  as  extensive  variables 

or  the  c - 1 values  of  their  ratios  as  intensive  variables.  For  two 
mixtures  with  the  same  gross  composition,  the  Q^'s  are  proportional  to 
each  other  and  the  proportionality  constant  is  the  ratio  of  the  masses 
of  the  two  systems.  Although  the  Q 's  are  the  natural  parameters  for 

K 

describing  the  gross  composition  of  the  mixture,  they  are  usually  unsuit- 
able for  computing  and  because  the  elements  are  present  in 

small  quantities  in  many  systems  of  interest.  It  is  therefore  convenient 
to  formulate  a general  method  for  describing  the  gross  composition  of  the 
mixture  so  that  the  conservation  of  mass  can  be  expressed  in  terms  of 


the  most  abundant  species  in  the  system. 

The  gross  composition  of  the  mixture  will  be  expressed  in  terras  of 
its  constituents,  and  the  constituents  sufficient  for  its  description 
will  be  called  the  components  of  the  system.  An  analytic  criterion  will 
be  formulated  for  selecting  a set  of  components.  The  criterion  is  based 
on  the  assumption  that  the  number  of  components  is  equal  to  the  number 
of  elements  c.  The  components  will  be  distinguished  with  an  index  j 
and  labelled  accordingly  with  j = 1,  ...,  c. 

Since  the  species  are  labelled  as  conntituents  with  i = 1,2,  . . .t . 

and  as  components  with  j = 1,2,  ...,  c,  the  index  i * i*(j)  is 

th 

introduced  to  show  that  the  species  labelled  as  the  j component  is 
also  labelled  as  the  i = constituent. 

It  is  convenient  to  rewrite  molecular  formulae  so  that  each  species 
in  the  system  can  be  represented  as  a vector.  The  molecular  formula  of 
the  constituent  labelled  i is  thus  written  as 
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where  X denotes  the  element  labelled  k,  and  a denotes  the  number 

lk 

of  atoms  of  element  k present  In  a molecule  of  species  i . A vector 
representation  of  each  species  is  then  readily  obtained  by  defining  the 
formula  vector  of  the  constituent  labelled  i as. 


y,  = (a.  , ...,ce  ) 
i ti  ic 


- 1.  ....  t 


(I-B-3) 


Consider,  for  example,  a system  containing  nitrogen,  oxygen,  carbon, 
and  hydrogen,  labelled  with  k * 1,  2,  3,  and  4 as  X1  = N , X3  = 0, 

X3  = C,  and  X*  = H.  Suppose  that  carbon  dioxide  is  labelled  with  i = 3 
and  water  with  i = 4.  Then  writing  Eq.  (I-B-2)  for  these  species  gives 


with 


and 


with 


C08  h y3  = NoOaCiHo 


Ha0  = Y = N?0xC3H8 


and  Eq . (I-B-3)  gives  their  corresponding  vector  representations  as 
y3  = (0, 2,1,0)  and  y4  = (0,1, 0,2). 

A necessary  and  sufficient  condition  for  the  selection  of  a proper 
set  of  components  is  that  the  formula  vectors  of  the  constituents  selected 
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as  components  be  linearly  independent.  This  condition  implies  that  the 
determinant  ^1  °*  formula  vectors  of  these  constituents  is 

nonzero.  Let  q denote  the  number  of  moles  of  the  j**1  component  in 
mass  Mo  of  the  hypothetical  system  consisting  of  components  only.  Then 
the  condition  for  the  hypothetical  system  of  components  to  have  the  same 
gross  composition  as  the  system  of  interest  follows  from  the  law  of 
conservation  of  masr.  as 


Z Of  q = o 

J=l  i*  (j)  k j 


k = 1 ,2 , . . . , c 


(I-B-4) 


The  quantities  q^(j  * ....  c)  are  called  the  stoichiometric  con- 

stants of  the  system  for  a particular  choice  of  components.  They  provide 
an  alternative  specification  of  the  gross  composition  of  the  system  to 
that  provided  by  the  Q^. 

A choice  of  components  is  generally  not  unique.  The  elements  always 
constitute  a proper  set  of  coaiponents.  And  in  the  case  that  they  are 
chosen  as  components,  Eq.  (I-B-4)  becomes 


E 5 q - Q 
J=i  Jk  j \ 


(I-B-5) 


where  is  the  Kronecker  delta  having  the  property  that  8^=1 

when  J = k , and  6 * 0 when  j ^ k,  However,  if  the  elements  are 

jK 

present  in  small  quantities,  as  is  usually  the  case,  it  is  more  convenient 
to  choose  components  from  the  more  abundant  species  in  the  system. 


We  are  now  in  a position  to  formulate  the  procedure  used  to  describe 
the  compos it tor  of  the  system  in  terms  of  components.  The  procedure  is 
based  on  the  fact  that  the  formula  vector  of  the  i constituent  can 
be  expressed  as  a linear  combination  of  the  formula  vectors  of  the 
components  as  follows; 


I -8-6 


<i  = 1,2,  , t) 


U-B-e> 


If  the  dissociated  elements  are  taken  to  be  the  components,  then  Eq.  (I-B-6) 

reduces  to  Eq.  (I-B-3)  and  » Cl£  The  expressions  in  Eq.  (I-B-6)  are 

i 3 i <3 

the  vector  representations  of  the  chemical  reactions  producing  the  con- 
stituents from  the  components;  the  corresponding  expressions  in  terms  of 
molecular  formulae  can  be  written  formally  as 


v1  = Ic  » I1*0' 

J=i  ^ij 


(i  = 1 ,2,  ...,t) 


(I-B-7) 


The  expressions  in  Eq.  { 1 -B-6)  contain  an  identity  for  each  constituent 

selected  as  a component.  Although  the  expressions  in  Eq.  (I-B-7)  can 

always  be  written  by  inspection,  using  the  customary  rules  for  balancing 

chemical  equations,  it  is  convenient  to  develop  a method  for  constructing 

the  matrix  that  can  be  used  in  the  computer.  The  coefficients  in 

the  g matrix  will  be  calculated  using  the  coefficients  a of  the 
Mij  ik 

formula  vectors  of  the  constituents. 

The  combination  of  Eqs.  (I-B-3)  and  (I-B-6)  leads  to  the  following 
matrix  equation 


<«ik>  * <*»,>< « 


ij  i*(j)k 


) 


(I-B-8) 


where  Of  „ . is  the  matrix  formed  from  the  formula  vectors  of  the 
i*(j)k 

components.  Since  is  by  definition  nonsingular,  Eq.  (I-B-8) 

can  be  inverted  to  give  the  equation 


%>  ‘ ‘V'V 


(I-B-9) 


^ “1 
where  (a  ) denotes  the  inverse  matrix  (a,  , ..)  . Equation  (I-B-9) 

kj  i*(j)k 
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can  be  rewritten  in  terns  of  the  coefficients  as 


$ 


ij 


A “iAj 


(O-b-io) 


Since  the  matrix  (0^)  ts  known,  Eq.  (I-B-9)  can  be  used  to  construct 
the  matrix  (g^)  a procedure  is  formulated  for  constructing  the 
matrix  (a^) * The  »atrix  (“ko  ) is  construe'  * by  selecting  an  appro- 
priate set  of  linearly  independent  formula  vectors  and  inverting  the 

corresponding  ( a , , , matrix.  The  Col  ) matrix  is  constructed  from 
i*(j)k  xj 

the  (a  ) matrix  with  a method  suitable  for  electronic  computers.  The 
ik 

selection  of  components  and  the  computation  of  the  (0^)  matrix  are 
carried  out  concurrently. 

The  rows  of  the  (a  ) matrix  are  added  to  an  initially  null  c x c 
ik 

matrix  one  at  a time.  After  each  such  addition  the  rows  are  tested  for 
linear  independence.  Linearly  independent  rows  define  components  and  * re 
retained,  linearly  dependent  rows  are  rejected,  and  the  process  is  continued 
until  a c x c matrix  is  obtained.  The  test  for  linear  independence  is 
made  by  beginning  the  reduction  of  the  matrix  to  triangular  form  (by  opera- 
tions on  rows)  with  the  first  such  addition  and  by  continuing  it  with  each 
subsequent  addition.  After  triangularizing , a row  is  linearly  independent 
if  it  contains  a nonzero  element  on  the  diagonal  or  to  the  right  of  the 
diagonal.  The  nonzero  element  is  placed  on  the  diagonal,  if  necessary, 
by  rearrangement  of  columns.  The  c x c matrix  is  used  to  construct 
the  ) matrix  by  back  substitution. 

The  (O  „ „ „ , (a  ) > and  (A  ) matrices  for  a system  containing 
i*(j)k)  ik  "ij 

C,  H,  O,  and  N (with  C,  COs,  Ha , and  Na  considered  as  components)  are 
given  below  to  exemplify  the  method  used  in  the  TIGER  code  to  describe 
the  composition  of  a thermodynamic  system; 
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ELEMENTS 


c 

H 

o 

N 

k = l 

2 

3 

4 

i*(J) 

C 

1 

1 

0 

0 

0 

co8 

3 

1 

0 

2 

0 

H* 

4 

0 

2 

0 

0 

n3 

6 

0 

0 

0 

2 

ELEMENTS 

c 

H 

0 

N 

(X 

ik 

k = 1 

2 

3 

4 

1 

C 

1 

1 

0 

0 

0 

CO 

2 

1 

0 

1 

0 

COz 

3 

1 

0 

2 

0 

Ha 

4 

0 

2 

0 

0 

HaO 

5 

0 

2 

1 

0 

Na 

6 

0 

0 

0 

2 

Oa 

7 

0 

0 

2 

0 

0 

8 

0 

0 

1 

0 

OH 

9 

0 

1 

1 

0 

H 

10 

0 

1 

0 

0 
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Section  I-B-3  presents  the  thermodynamic  description  of  nonideal 
heterogeneous  systems  used  in  the  TlfTER  code.  Systems  in  mechanical  and 
thermal  equilibrium,  with  condensed  phases  considered  as  pure  substances 
to  exclude  solid  and  liquid  solutions,  are  treated  under  the  assumption 
that  an  equation  of  state  exists  for  each  phase.  The  equations  of  state 
are  restricted  only  by  thermodynamic  identities  and  stability  conditions. 
This  description  is  convenient  to  use  in  calculations  of  the  thermodynamic 
state  of  such  systems  at  one  of  the  following  points,  tp,T)  , (p,h) , (p,s) , 
(v,T) , (v,e) , (v,p)  or  (v,s) , when  the  composition  is  in  chemical  equilib- 
rium, in  partial  equilibrium,  or  completely  frozen  in  metastable  equi- 
librium. In  the  treatment  of  partial  equilibrium,  part  of  the  composition 
is  prescribed  to  be  frozen  and  the  remainder  is  assumed  to  be  in  chemical 
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equilibrium.  Calculation  'A  the  state  at  a point  Involves  calculating 
the  composition  and  evaluating  the  other  state  variables  with  thermo- 
dynamic identities.  As  stated  in  Section  l-B-2,  the  composition  is 
restricted  a priori  to  s gaseous  species  with  mole  numbers 
n^(i  - 1,2,  ...,  s),  t - s condensed  species  with  mole  numbers 
#*  (i  s 5 + 1,  . . . , t)  , and  oust  satisfy  the  stoichiometric  conditions 
expressing  the  law  of  conservation  of  mass  for  th?  system. 

Calculation  of  the  state  is  easiest  when  the  composition  if  completely 
frozen  because  the  mole  numbers  can  be  chosen  to  satisfy  the  stoichiometric 
conditions.  The  calculations  are  more  complicated  for  the  equilibrium 
cases  because  the  composition  must  satisfy  the  equilibrium  conditions  as 
well  as  the  stoichiometric  conditions.  Since  calculation  of  the  equilib- 
rium composition  is  based  a priori  on  the  assumption  that  some  of  the  con- 
densed species  are  present  and  the  remainder  are  absent , the  results  of 
the  calculation  must  be  used  to  test  if  this  assumption  satisfies  the 
equilibrium  conditions  for  condensed  species.  To  satisfy  these  equi- 
librium conditions,  the  mole  number  of  a species  i assumed  to  be 

% 

present  must  satisfy  the  condition  n^  s 0,  and  the  chemical  potential 
of  a species  assumed  to  be  absent  must  exceed  its  chemical  potential  in 
the  gaseous  phase. 

a . The  Stoichiometric  Conditions 

The  stoichiometric  conditions  expressing  the  conservation  of 
mass  for  the  system  will  be  formulated  by  expressing  the  formula  vectors 
of  the  constituents  in  terms  of  the  formular  vectors  of  the  components. 

A set  of  components  is  first  chosen,  and  the  corresponding  stoichiometric 
constants  are  evaluated.  The  system  is  th  n transformed  into  a system 
containing  only  these  components  by  expressing  the  conslitutents  as 
linear  combinations  of  them.  T!v;  stoichiometric  conditions  are  then 
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obtained  by  equating  the  mole  numbers  of  the  components  in  the  transformed 
system  to  tbeir  stoichiometric  constants.  Following  this  procedure  with 
the  notation  introduced  previously  in  Section  l-B-2  gives  the  stoichio- 
metric conditions  as 


A Vi + A-n  V*  * s J * 


(I-B-ll) 


It  is  convenient  to  introduce  the  index  i '(ra>  , with  m - 1,2,  ....  t - s, 

so  that  the  condensed  species  labeled  as  m is  also  labeled  as  the 
th 

i = i 'tm)  constituent.  It  is  also  convenient  to  denote  the  mole  numbers 

of  the  condensed  species  by  N so  that 

81 


i'(w) 


BS  - 1,2, 


(I-fi-12) 


The  condensed  constituents  assumed  to  be  present  and  in  equilibrium  with 

the  gaseous  phase  will  always  be  Chosen  as  components  as  in  the  original 

veision  of  the  TIGER  code  formulated  by  R.  S.  Brinkley,  Jr.  It  is 

necessary  to  set  i*(j)  * i '(sn)  , and  n*/  . = N for  ra  = j = 1,2 p 

H«)  J 

to  make  the  i'  map  compatible  with  the  i*  map  introduced  previously 

for  mapping  components  into  constituents.  The  condensed  species  assumed 

to  be  frozen  and  not  in  chemical  equilibrium  with  the  gaseous  phase  will 

be  labeled  with  m ~ p + 1,  ....  p',  and  the  remainder  that  are  assumed 

to  be  absent  (but  to  satisfy  the  equilibrium  conditions)  will  be  labeled 

with  m=p'+l,  .«.,  t — s . Thus  the  mole  numbers  of  the  condensed 

species  are  subject  to  the  following  restrictions:  N = N ^0  for 

m j 

a 5 J ® 1,2,  ...,  p;  H £ 0 for  m = p + 1,  ....  p';  and  N =0  for 

m n 

a * p'  + 1,  ....  t - s.  The  summation  term  over  the  condensed  species 
In  iSq.  (I-B-ll)  can  therefore  be  rewritten  as 


£ fln*=?dN  + 2 A v N 

i-s+i  HlJ  i m—i  mj  m m=p+i  *“i  (m)  j in 


3 = 1 C (I-B-13) 


It  is  convenient  for  comparative  purposes  to  wtii.*?  an  equation 
similar  to  (I-B-13)  for  the  gaseous  summation  term  in  Eq.  (I-B-il) , 
although  the  gaseous  species  will  not  be  treated  with  this  equation  in 
the  code.  Assume  that  s'  of  the  gaseous  species  are  in  equilibrium 
while  the  remaining  s - s'  are  frozen,  and  introduce  the  index  f(g) 
with  g = 1,  b to  indicate  that  the  gaseous  species  labeled  as  g 

is  also  labeled  as  the  i - i(g )tb  constituent.  The  gaseous  species 
assumed  to  be  components  will  be  labeled  with  g - 1,  ...»  c - p,  those 

assumed  to  be  in  equilibrium  with  g = c - p + 1,  and  those 

assumed  to  be  frozen  with  g * s'  + 1*  As  in  the  treatment  of 

the  condensed  species,  it  is  necessary  to  set  i{g>  = i*(g  + p)  for 

g = 1,  ...»  c - p to  make  the  i map  compatible  with  the  i*  map  for 
identifying  components  as  constituents.  With  this  notation  the  summation 
term  over  the  gaseous  species  in  Eq.  (I-B-ll)  can  be  rewritten  as 


8.  n. 
Hi  J i 


r 6 n»  . 

g=i  g+p,j  i(g) 


s'  s 

gSc-pt-i®i(g)  jn{(g)+  g=s'+1^f(K>  jnf<e> 


0=1. 


(I-B-14) 


Since  m - 1,  ...»  p and  ■ = p'  + 1,  ....  t-s  are  used  to  label  the 
condensed  species  assumed  to  be  in  chemical  equilibrium,  and  m = p + I,  pJ 

are  used  to  label  the  condensed  species  assumed  to  be  frozen,  the  condition 
for  all  the  condensed  species  to  be  in  equilibrium  is  p'  - p,  and  the 
conditions  for  them  all  to  be  frozen  are  p = 0 and  p'  = t - s.  Similarly, 
since  g = 1,  ...,  s'  are  used  to  label  the  gaseous  species  in  chemical 
equilibrium  and  g = s'  + 1,  ...,s  are  used  to  label  the  frozen  species, 
the  condition  for  all  the  gaseous  species  to  be  in  equilibrium  is  s = s ' , 
and  the  condition  for  them  all  to  be  frozen  is  s'  = 0.  Thus  p = p'  and 
s = s'  when  the  composition  of  the  heterogeneous  system  is  in  complete 
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f 


i--?a:^sp¥w-.'U4  rv?*atsr*<’ 
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equilibrium,  and  p = s * =0  and  p*  = t - s when  it  is  completely 
frozen.  When  t = s,  no  condensed  species  are  present  and  tlic  system  is 
a homogeneous  gas  phase.  Furthermore,  when  t = s = s',  the  composition 
of  the  gaseous  mixture  is  in  chemical  equilibrium,  and  when  t - s = s'  ~ 0, 
it  is  completely  frozen. 

In  calculations  of  the  thermodynamic  state,  the  mole  numbers  of 
the  equilibrium  species  are  calculated  using  the  stoichiometric  conditions 
and  the  equilibrium  conditions.  The  c - p gaseous  components  are  chosen 
first  from  the  gaseous  constituents  assumed  to  be  in  equilibrium;  they  must 
include  frozen  species  when  s*  < c - p,  and  must  be  chosen  from  the 
frozen  species  when  s/  = 0.  It  is  convenient  for  computational  purposes 
to  combine  Eqs.  (I-B-ll)  and  (I-B-13)  and  rewrite  the  stoichiometric 
conditions  as 


m=p+i ^i'(m)  1N  m i=i ^i 


+ N = R, 


J =1,2, 


(I-B-15) 


E'  s 

£ 9/r>N  + r;*n=q 

m=p+iai  (m)j  m i=i°ij  i j 


j * p + 1,  . ..,  c (I-B-16) 


These  equations  are  generalizations  of  the  stoichiometric  conditions 
formulated  by  Brinkley  for  the  case  of  equilibrium,  and  thus  they  reduce 
to  the  stoichiometric  conditions  presented  in  the  original  documentation 
of  the  TIGER  code  when  p = p7  and  the  first  summation  term  vanishes. 

For  the  homogeneous  system  when  p/  = p = 0,  Eq.  (I-B-15)  vanishes  and 
Eq,  (I-C-16)  reduces  to  the  stoichiometric  conditions  for  the  gaseous 
phase . 

Let  n denote  the  total  number  of  moles  of  gas  in  our  mixture 
of  mass  Mo  so  that 

n = I n O-B-17) 

i=l  i 


I-B-14 


It  is  useful  for  computational  purposes  to  have  another  expression  for 
n and  such  a relationship  can  be  obtained  by  combining  Eqs.  iI-B-16) 
and  (I-B-17).  Thus,  summing  Eq.  (I-B-16)  over  j and  subtracting  the 


from  Eq.  (I-B-17)  leads  to  the 

equation 

s 

P' 

n = q + Z (1  - a )n.  - 
g t=i  i 

m=p+i^i  '<m)Nm 

(I-B-18) 

qg’  * and  ^i(m)  are  deflned  b?  the  equations, 

c 

q = £ q 

S J*P+l  J 

c 

(I-B-19) 

J'p+l  ®ij 
c 

(I-B-20) 

'i'On)  = j^B+l  8i'(m)j 

(I-B-21) 

It  is  also  useful  to  define 
: in  terms  of  the  components  as 

c 

the  reference  mass  Mq 

of  the 

Mq  = £ q M 

J=i  J 3 

(I-B-22) 

t h 

where  is  the  molecular  mass  of  the  j component,  and  to  normalize 
the  system  by  choosing  Mo  as  100  grams.  Then  the  mass  balance  equation 
for  the  normalized  system  is 


100  * £ (100  q /M0)M 

J-*  J J 


(I-B-23) 


and  the  corresponding  stoichiometric  coefficients  are 


q ' = 100  q /Me 

3 3 


(I-B-24) 
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Thus,  the  equations  for  the  normalized  system  are  obtained  from  those  for 

the  reference  system  of  mass  Mo  by  replacing  q by  100  q /Mq . 

>3  3 

If  we  consider  tl»e  decomposition  products  of  TNT,  C^fcO^N.,  , 
as  an  example  and  choose  C,  CO,  Ha , and  N-,  as  components  according 
to  the  reaction 


C^HgOgNg  _ C + 6 CO  + 5/2  H2  + 3/2  Na 


then  the  normalized  stoichiometric  coefficients  are  q'  = 100/M,  qr  - 

c CO 

600 Ai,  q^  = 250/M,  and  ^ 150/M,  where  M = 228.18  is  the 

molecular  weight  of  TNT. 

b.  Equations  of  State 

The  gas  phase  is  treated  with  an  equation  of  state  of  the  form. 


p = ptp.T.ni ng) 


(I-B-25) 


expressing  the  pressure  as  an  explicit  function  of  the  gas  density  p, 
the  temperature,  and  the  mole  numbers  of  the  gaseous  constituents.  It  is 
convenient  to  introduce  the  variable  p through  the  identity 


fi-  = fi- 

Mq  M 


(I  -B-26) 


so  that  l/$  is  the  volume  of  gas  per  unit  mass  of  mixture,  and  p = p 

when  M = Mo  and  no  condensed  species  are  present.  Equation  (I-B-25) 
8 

is  then  written  more  explicitly  as 


p = , ...,  n) 

Mo  s 


(I-B-27) 


with  $ an  imperfection  term  that  approaches  1 as  the  gaseous  mixture 


becomes  ideal. 
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„ th 

The  chemical  potential  of  the  i gaseous  constituent 


is  written  as 

j^/RT  - p°/*T  + ^ + 


c+i 


i = X,  ... , s (I -8-28) 


where  - ja® <T>  is  the  chemical  potential  of  the  i~"  gaseous  con- 
stituent in  its  standard  state  at  unit  pressure,  and  , r ^ , and 


th 


r are  defined  by  the  equations 


«i  * J""i 


(I  -B-29) 


r - In  (RTp/Mo) 
c+i 


(l-B-30) 


A 


r * 

i Jq 


PK&-V  i % 

Un&bn)  Jfi 


(i-B-31) 


with  the  integration  performed  along  an  isotherm.  The  imperfection  term 
T can  be  considered  as  the  logarithm  of  an  activity  coefficient  and  can 
be  written  in  terms  of  ♦ as 


* 

A 


_ f*"  - . a*  ,1  an 

r - i . ♦ + n-  - 1 

i *q  t.  3t>.  J A 
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The  hermodynamic  identity 


d Jtn  T \RT  J 


x 

RT 


i = 1, 


(I-B-33) 


is  used  to  introduce  the  molar  enthalpy  IT  and  the  reduced  molar  enthalpy 
of  the  ith  gaseous  constituent  in  its  standard  state. 


The  condensed  species  are  treated  with  equations  of  state  of 


the  form 
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II 

CT,p> 

i = s 

+ 1 

t 

(I-B-34a) 

(T.p) 

( I -B-34b) 

* it 

where  and 

denote  the  molar 

volume 

and  molar 

enthalpy  of  the 

.th 

i condensed  constituent.  It  is  convenient 

to  label  these  molar  quanti- 

ties  with  the  index 

m introduced  earlier  for  condensed  species 

, so  that 

**i'«  = 

V* 

m 

m = 1 

t • • • > f “ 

s 

(1 -B-35a) 

“t'w  * 

H* 

m 

(I  -B~35b) 

and  i'(ra)  = i*(j) 

for  m » j « 1 ,2 , 

• • . > p. 

Use  will 

be  made 

of  the 

derivatives 

,b  **•  v«  X 

Of  (T,p) 
m 

- ( n ) 
\B  fn  T / 

P 

(I-B-36a) 

£ <?.p> 

a in  V* 

■ {rr.T-\ 

(I-B-36b) 

and 

dH* 

t * \ 

V srr  ) 

p 

c*«. 

(I-B-37a) 

II 

jfljl 

v*  a - a*) 

m m 

( I -B  -37b) 

The  derivative  c?  is  related  to  the  coefficient  of  thermal  expansion 

HI 

* 

at  constant  pressure  and  a to  the  coefficient  of  compression  at 

th  . * 

constant  temperature  of  the  m condensed  species;  C is  the 

pm 

corresponding  constant  pressure  molar  heat  capacity. 

The  chemical  potentials  of  the  condensed  species,  in  contrast 
to  the  chemical  potentials  of  the  gaseous  species,  are  functions  only  of 
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temperature  and  pressure  and  will  be  written  as  either 


Pt  ” »P> 


i = s + 1 , . . . , t 


<i  -»-;$«) 


P„  - H*<T,p> 
a a 


ffl  ”■  1 f « * a i t “ S 


(I  -B-39) 


with  i'(ra)  = i*(j)  and  y . , * p.*  for  m = j * 1,  ...,  p.  The 

i (tt)  m 

cnemical  potentials  of  these  species  in  the  gaseous  phase  are  denoted  by 

with  i = s + 1,  ...»  t.  It  is  important  to  note  that  the  index  i 

on  the  gaseous  potential  has  the  range  i = 1,  ...»  t and  covers 

all  the  constituents,  even  though  the  species  labeled  with  i = s + 1,  t 

are  assumed  a priori  to  be  condensed  constituents.  It  is  convenient  for 

computational  purposes  to  define  the  reduced  molar  volume  ® and  the 

m 

reduced  molir  enthalpy  ^ of  the  condensed  species  with  the  thermodynamic 
identities 


P*(P,T)\  pV 


&p\  RT 


m _ * 

RT  " m 


a K(P'T,I  “■ 


1 dTl  RT 


(I-B-40) 


(I-B-41) 


Equations  (I-B-26)  , (I-B-40),  and  (I-B-41)  are  used  to  obtain 
a convenient  expression  for  the  specific  volume  of  the  mixture.  The 
volume  V of  the  mixture  is  written  as  the  sum  of  the  volumes  of  the 
phases  as 


V=Mo/A+ZNV*+  ? N V* 
J=i  J j m=p+i  m m 


(I-B-42) 


by  remembering  that  m = j for  m = 1,  ...,  p.  The  equation  for  the 
specific  volume  is  obtained  as 


y = 1/p  + JL  { | N (ff  «■  ? N <1?) 

Ro P ' J=i  J J m=p+i  m m J 

by  dividing  Eq.  (I-B-42)  by  Mq . 


( I -B-43) 
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4. 


Calculation  of  the  Equilibrium  Composition 


a . Equilibrimi  Conditions  and  Iteration  Parameters 

It  1b  convenient  to  consider  the  case  of  complete  equilibrium 
when  p = p',  s = s',  and  the  mixture  contains  no  frozen  constituents. 
Then  equilibt  ' ura  conditions  must  be  formulated  for  the  s gaseous  consti- 
tuents, for  the  p condensed  constituents  assumed  to  be  present,  and  for 
the  remaining  t - <s  + p)  condensed  constituents  assumed  to  be  absent. 
Use  will  be  made  of  the  fact  that  the  chemical  potentials  of  the  species 
in  equilibrium  satisfy  the  same  equations  as  their  formula  vectors. 
Equation  (I-B-6)  is  thus  used  to  express  the  formula  vectors  of  the  con- 
stituents as  linear  combinations  of  the  formula  vectors  of  the  components. 
The  equations  for  the  formula  vectors  are  written  as 


y 


1 


-P  + rc 

J=i  *ijyi*(j>  j=p+l^ijyi*( j) 


i = 1 , ....  t (I-B-44) 


because  the  p condensed  constituents  assumed  to  be  present  must  be  chosen 
as  components.  Replacing  y by  u in  Eq.  (I-B-44)  and  taking 

account  of  Eqs.  (I-B-38)  and  (I-B-39)  leads  to  the  equation 


_ * e V. 

Hi  + j=p+i^i J^i*( j)  i =1,  ....  t (I-B-45) 

Equation  (I-B-45)  expresses  the  chemical  potential  of  < ach  constituent  in 
the  gaseous  phase  in  terms  of  the  chemical  potentials  of  the  condensed 
and  gaseous  components.  Thus  for  i = 1,  ....  s it  is  an  expression  for 
the  gaseous  constituents,  but  for  i ■ s + 1,  ...,  t it  is  an  expression 
for  the  gaseous  species  assumed  a priori  to  be  condensed  constituents. 

Equation  (I-B-45)  contains  c identities  for  the  values  of 
i * i*(j)  labelling  the  constituents  chosen  as  components.  Introducing 
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jy  = 1,  ...»  c as  a dummy  component  index  so  that  /}  * 6 , , and 

i*(j  )j  J J 

setting  i = i*(j'}  in  Eq.  (I-B-45)  gives  these  identities  as 


“i  = V(j) 


J = 1*  P 


(I -fi-46a) 


Ai*(j) 


3 = p + 1 c (I  -B-46b) 


Whereas  Eq.  (I-B-46a)  is  the  condition  for  a condensed  constituent  that 
is  present  to  be  in  equilibrium  with  the  gaseous  phase,  Eq,  (I~B~46b)  is 
an  identity  specifying  the  gaseous  constituents  chosen  as  components.  We 
must  now  consider  the  condensed  constituents  assumed  to  be  absent  with 
i C s + 1,  ...s  t and  i = i*(m) , m = p+  1,  . t - s.  The  equilibrium 
condition  for  each  of  these  species  is  that  its  chemical  potential  in  the 
condensed  phase  must  be  greater  than  its  chemical  potential  in  the  gaseous 
phase.  Gibbs'  conditions  for  chemical  equilibrium  in  the  heterogeneous 
mixture  can  thus  be  written  as 


Ai  j=l^ij^j  + j=p+l^ijfli*(j) 


i = 1, 


(I-B-47a) 


**i*Cj) 


J = 1,  p 


(I  -B-47b) 


^i'(m)  j=i^i  '(ra)  + j=p+i^i '(m)  jHi*<j) 


m — p + 1 , . . . , t 


Thus  the  gaseous  constituents  in  an  equilibrium  mixture  must  satisfy  Eq. 
(I-B-47a) , the  condensed  constituents  assumed  a priori  to  be  present  must 
satisfy  Eq.  (I-B-47b) , and  those  assumed  to  be  absent  must  satisfy 
Eq.  (I-B-47C). 
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i 

1 


When  * complete  equation  of  state  of  the  mixture  is  known, 

Eqs.  (IHS-47)  and  the  stoichiometric  conditions  obtained  by  setting  p = 
p * in  Eqs.  (I-B-15)  and  (IHB-16)  are  sufficient  for  calculating  the 
equilibrium  composition  at  any  one  of  the  following  specified  points 
(p,  T),  (p,  h)  , (p,  s>  , <v,  T)  , (v,  e)  , <v,  p)  , or  (v,  s)  . Since  ^ 
is  known  as  a function  of  the  n^'s,  the  s - c + p independent  equilib- 
rium conditions  together  with  the  c stoichiometric  conditions  constitute 
a set  of  s + p equations  for  calculating  the  s + p unknown  mole  numbers 

n (1  + 1,  s)  and  N ( J = 1 , ....  p) . The  computed  compositions  must 

— 3 

be  checked,  however,  to  test  the  validity  of  the  a priori  assumptions  that 

had  to  be  made  about  the  presence  (K  & 0 for  j * 1 , . . . , p)  and  the 

3 

absence  (N  = 0 for  a = p + l,  ...,t-s)  of  condensed  constituents. 

m 

The  assumptions  are  valid  when  the  calculated  composition  satisfies  the 
condition  a 0 and  Eq.  (I-B-47c),  and  they  are  invalid  when  the  calcu- 
lated composition  does  not  satisfy  these  conditions.  In  the  latter  case, 
the  assumptions  must  be  changed  and  the  composition  must  be  recalculated 

until  N a 0 and  Eq.  (I-B-47c)  is  satisfied. 

3 

The  equations  used  in  the  TIGER  code  to  compute  the  equilibrium 

composition  will  now  be  derived.  The  first  step  is  to  obtain  an  expression 

for  the  equilibrium  conditions  of  the  gaseous  species  in  terms  of  their 

mole  numbers.  Thus  substituting  u and  u given  by  Eq.  (l-B-28) 

* 

into  Eq.  (I-B-47a)  and  making  use  of  Eq.  (I-*.  ?9>  leads  to  the  following 


equation  for  E.  - tn  n (i  = 1,  . . . , s) , 
i i 


4 = - g - r + £ n*/RT  + I i n,  + (fi4  - 1)1 

Hi  i i J=rlj  j J=p+1  i J J 


(I-B-48) 


where 


Si  = *1™  ~ 


( I -B  -49) 


r = e + r 

j n*<j)  i*  ( j) 


J = p + 1 c 


(1 -B-5G) 
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Equation  (IHB-48)  reduces  to  Eq.  (I-B-50)  for  the  gaseous  constituents 


l 

£ 


chosen  as  components  because  A , ,,  = 6 , g - O , and  A = 1 for 

i*<3  >)  JJ  i pi 

i = i’s'(j/)(j/  = p + 1,  . . . , c)  . The  next  step  is  to  express  the  inequali- 
ties of  Eq.  <i-B~4?c)  in  terms  of  the  mole  numbers  of  the  gaseous  species. 
The  combination  of  Eq.  (I-B-28)  end  (I-B-47c)  leads  to  the  expression 


*m  3=  i^i'tmJj^j  j=p+i^i '(m)  3 j + '(ra)  c+i  < ° 


where 

* 

gm 


for  » - p r 1 , 


* 

u.  /RT  - 
m 


t - s 


C I -B  -51) 


( I -B  -52) 


It  is  now  convenient  to  simplify  Eqs.  (I-B-48)  nod  (I-B-51)  by  introducing 
the  parameters 
* 

+1  = p /RT  J = 1,  ....  p (I-B-53) 

J J 

and 

/»!,.*  ’ *1  - 1 


For  notational  convenience,  the  index  j used  for  the  - parameters  will 
be  extended  from  c to  c + 2 so  that  Eq.  (I-B-48)  can  be  written  as 

«i  ■ - «1  - ri  ♦ "i  »iy3  5 <<‘B‘55) 

and  Eq.  (I-B-51)  can  be  written  as 
c+i 

* 

- g + " A.  ,,  xr+r  <0  m = p+  l,...,t-s  (I  -B-56) 

m j=i  (nOj  j c+i 

Other  equations  needed  to  compute  the  equilibrium  composition  are  the 
stoichiometric  conditions 


I 

j 


I 


I 

I 
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' Wr  . J ■>,.  -He.:  ;» 


i 

i 

\ 

; 

* 

5 

t 

l 

f 

f. 

t 


f.  A n + N = q 

i=x°lj  i J 0 


J~l»  • • • i p 


( I -B-57) 


$ 

x.  « n = q 


J = p + 1 , . . . , c 


(I-B-56) 


obtained  by  setting  p = p'  in  Eqs.  Cl-B-15)  and  (I-B-16). 

Equations  (I-B-53),  (I-B-55),  (I-B-57),  and  (I-B-58)  can  be  used 
to  compute  the  equilibrium  composition  with  a set  of  condensed  components 
specified  a priori  by  the  assumption  about  the  presence  of  condensed  con- 
stituents, and  then  Eq.  (I-B-56)  can  be  used  to  test  whether  this  assumption 

is  valid.  It  is  convenient  to  define  r>  by  the  equation  I 

C+3  j 

I 

! 

t,m.s  “ tnT  a -*-59>  1 


For  fixed  values  of  and  ^ , Eqs.  (I-B-53),  (I-B-55),  (I-B-57)  , 

and  (I-B-58)  constitute  (c  + p + s)  equations  for  the  (c  + p + s) 

unknown  quantities  ~ (j  = 1 , . . . , c) , N (j  » 1 , . . , p) , and  n. (i  - 1, 

J J i 

...,  s)  . Equations  (I-B-53),  (I-B—  55)  , and  (I-B-58)  can  be  used  to 

determine  the  c + s values  of  r and  u. , however,  because  K 

J ^ J 

appears  only  in  Eq.  (I-B-57) . Solving  the  former  equations  thus  determines 

the  composition  of  the  gas  phase  at  fixed  temperature  and  density, 

provided  that  the  condensed  species  assumed  to  be  present  satisfy 

the  inequalities  H a ° (j  = 1,  p) . And  solving  the  latter  equation 

for  the  N 's  with  the  calculated  gaseous  composition  determines  the 
J 

composition  of  the  condensed  phases.  According  to  Brinkley,  this  method 
of  treating  the  equilibrium  problem  in  terms  of  gaseous  and  condensed 
components  is  consistent  with  the  Gibbs  Phase  Rule, 

An  iterative  procedure  is  required  for  the  simultaneous  solution  of 
the  equilibrium  and  stoichiometric  conditions  because  Eq.  (I-B-58) 
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is  linear  in  the  mole  numbers  of  the  gaseous  species,  but  Eq.  (I-B-55)  is 
transcendental  in  these  mole  numbers  and  their  logarithms.  The  n parameters 
are  considered  as  independent  variables  in  Brinkley's  iterative  procedure 
for  solving  the  equilibrium  problem  at  a specified  point.  It  should  be 
remembered  that  the  index  j originally  introduced  for  components  has  the 
range  1=1,  . c + 2 when  used  for  the  *•  parameters.  That  the 
*■  {J  =1,  ....  c + 2)  parameters  can  be  chosen  as  independent  variables 
follows  from  the  equations  presented  earlier  in  this  volume.  Specifically, 

Eq.  (l-B-55)  can  be  rewritten  as 


G.  ( 

l 


r ) 

C+2 


with 


- 8, 


c+i 


+ 8.  ~ 

u J 


( I -B  -60) 


(IHB-61) 


because 
and  inT 


defined  by  Eq . (I-B-49)  is  a function  of  temperature  only, 

~ ' Equation  (IHB-60)  can  then  be  regarded  as  an  expression 

= n^,  ....  'c+a)  (I-B-62) 


defining  the  mole  numbers  of  the  gaseous  species  as  implicit  functions  of 

fhe  ” because  - in  and  is  an  explicit  function  of  p,  T, 

and  n , 
i 

r = r(nx,  ....  ns.  rc+1.  rc+a)  (I-B-63) 


It  follows  from  Eq,  (I-B-62)  that  Eq.  (I-B-57)  can  be  regarded  as  an 
expression 


- ) 
C+E 


(I-B-64) 
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i 


defining  the  mole  numbers  of  the  condensed  species  as  implicit  functions 

of  the  ■’r  parameters. 

J 

Brinkley's  approach  to  the  equilibrium  problem  was  to  construct 
a set  of  c + 2 equations 

F (■%,...,  r>  > - 0 j = 1,  .. . , c + 2 (I-B-65) 

J 1 C+2 


so  that  the  set  of  satisfying  these  equations  defines  the  equilibrium 

3 

compos i t ion . This  set  of  ^ that  satisfies  Eqs.  (I-B-65)  will  be  called 

the  set  of  equilibrium  values  of  ” . It  is  convenient  to  write  the  itera- 

3 

tive  equations  used  to  solve  (I-B-65)  for  the  equilibrium  values  of  ^ 
with  the  Kewton-Raphson  technique  as 


^ / l 

k?i  ' S'  V'1  ’ 

k 


c+2 


(I-B-66) 


i+i  i i t'"> 

where  - *•_  - ”,  and  +>  is  the  i ac  roxiraation  to  the  equi- 

k k k k 

librium  value  of  *•  . It  should  be  noted  that  the  index  k introduced 
k 

originally  for  atoms  is  also  being  used  as  a dummy  index  for  the  iteration 
parameters  in  order  to  simplify  the  notation  in  the  remainder  of  the 
documentation.  The  partial  derivatives  &F  must  be  known  as  explicit 

functions  of  the  parameters  in  order  to  generate  a solution  to  Eqs. 

(I-B-65)  with  (X-B-66) . This  being  the  case,  an  initial  approximation  to 
the  solution  (t-1  , . . . , ^ ) is  used  to  generate  a second  approximate 

1 C+2 

3 3 

solution  (T\  > ...»  r „)  with  (I-B-66),  and  the  procedure  is  continued 
TL  c+* 

until  preassigned  conditions  for  convergence  to  the  solution  are  satisfied. 

The  conditions  used  to  test  for  convergence  in  the  present  version  of  the 
P 

code  are  Z Fa/p  < Cx  , I F*7(c  - p>  < (a,  Fa  < <Tt  , and  F2  , < f4 
3=1  J 3=P+l  J c+i  c+2 

where  f x , f s . (3  > and  are  preassigned  small  numbers.  The  Eqs. 

(I-B-53)  and  (1-8-58)  were  chosen  to  define  c of  the  expressions  in 
( I — 8 -65)  as 
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F » p /RT  - - = 0 J * i,  ....  P (I-B-67) 

J J 3 


s 

F - <*  -Iji.  = o j = p + l,  c (i-u-6h) 

0 3 1=1  1 


and  the  remaining  two  equations 


F 

e+i 


F 

c+a 


0 

0 


(I-B-69) 


were  used  to  define  the  constraints  imposed  by  specifying  the  two  indepen- 
dent state  variables  at  the  thermodynamic  point  of  interest.  The  differ- 
ent expressions  for  F and  V for  the  seven  different  point 

c+i  c+a 

calculations  considered  earlier  will  be  presented  later.  A double 

iterative  scheme  is  required  to  solve  Eqs.  (I-B-67)  and  (I-B-68)  for  the 

equilibrium  values  of  + because  Eq.  (I-B-60)  must  be  solved  for  n 

in  order  to  evaluate  the  F *s  and  the  partial  derivatives 

3 


s ^ 
i=l^ij  ni 


J * p + 1,  ...,  C (I-B-70) 


obtained  from  the  stoichiometric  conditions  expressed  by  (I-B-68). 

The  equilibrium  values  of  - are  determined  in  an  outer 

J 

iteration  x.ith  Eqs.  (I-B-67)  to  (I-B-69),  using  values  of  n calculated 
in  an  inner  iteration  with  Eq,  (I-B-60)  and  the  Newton-Raphson  technique. 
The  equations  in  (I-B-60)  were  rewritten  as 


- r,  * V. 


71  ) - 0 

c+a 


(I-B-71) 


i — 1 , ...»  s 

and  solved  for  ^ = jtn  to  determine  the  mole  numbers  corresponding 
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I 

3f 


to  a fixed  set  of  the  ^ parameters  in  the  inner  iteration.  The  itera- 
tive equations  used  to  solve  (I-B-71)  were  written  similarly  to  those  used 
to  solve  (I-B-65)  as 


s 

T— 

r=i 


with 

aw . ar 

i = _ g _ ( — i n 

bP  ir  J r 

r 


ifk  = ' Mq  k“p  \ dp 

dn^  J0RTpi5n  Sn  Ip 


i - 1,  . . , , s (I-B-72) 


(l-B-7'i) 


k+l 

r 


and  £ the 

’ v* 


th 


approximation  to  P 


When  the  equilibrium  composition  of  the  gaseous  phase  has  been 
determined,  the  mole  numbers  of  the  condensed  phases  are  calculated  with 
Eq . ( I -B -57) , and  Eq.  (I-B-56)  is  used  to  test  the  validity  of  the  a priori 
assumption  about  the  condensed  phases  that  is  required  to  perform  the 
equilibrium  calculation.  When  the  equilibrium  composition  of  the  mixture 
is  known,  the  thermodynamic  state  at  a specified  point  is  computed  with 
thermodynamic  identities. 

b.  State  Equations 

The  equations  will  now  be  presented  for  computing  the  thermo- 
dynamic state  in  the  seven  different  point  calculations  after  the  equi- 
librium composition  has  been  determined.  The  pressure  p , the  specific 
volume  v , the  specific  enthalpy  h , the  specific  Gibbs  free  energy 
g , the  specific  internal  energy  e , the  specific  entropy  s , and 
their  derivatives 


I-B-28 


will  all  be  regarded  as  Implicit  functions  of  the  iteration  parameters 


5 

t 

* 

;• 


f-. 


. , r,  ) . The  partial  derivatives  of  these  variables  with  respect 
c+2 

that  are  used  in  the  iterations  will  be  presented  later. 


The  equation  for  the  pressure  is  written  as 

A 


P 


. 252£  *r  , 

Mo  1 


-r  ) 
C+2 


(I-B-74) 


with  0 regarded  as  an  implicit  function  of  the  iteration  parameters 

instead  of  an  explicit  function  <J>  = $(£,  T,  n,  , . ..,  n ) specified  by 

i s 

the  equation  of  state  of  the  gaseous  mixture.  The  equation  for  the 
specific  volume  of  the  equilibrium  mixture  is  obtained  by  setting  p = p# 
in  Eq.  (IHB-42)  as 

v(-l  , . ..,  r,  ) = 1/p  + 1/Mo  ( £ NV*)  (I-B-75) 

1 'c+2  J=1  J J 


The  specific  enthalpy  h of  the  mixture  is  given  by  the  expression, 

P * s A 

M0h  = IS  H + .Z  n r + pMo/p  + RTCc  - n)  (I-B-76) 

0 1 0 0 l\ii 


with  the  imperfection  term  < defined  by  the  following  integral  along 
an  isotherm 


(I-fi-77) 


Another  expression  for  ht^ , 
is  obtained  as 


T]  ) with 
'C+2 


ft  = v°  - r ft  v0 

°i,c+3  *i  j=p+iaijAi*(j) 


Moh  = N H 
0 3 


+ T.  q H° 

0=p+i  j i*(j) 


+ RT 


I fl  t 
i=lMi ,c+a 


A 

P«o  /P 


+ RT(C  - n)  (I-B-78) 


by  multiplying  Eq.  (I-B-58)  by  and  adding  the  resulting  expression 

i*  C J) 

to  Eq.  (I-B-76).  The  frozen  heat  capacity  of  the  mixture  is  given  by 


I -B  -29 


the  equation 


P * c 

c = IK  + EqC-  , v 
f J=*  3 PJ  J=p-^l  J Pi*(j) 


s 

v 
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(£  C°)n 
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n) 


(I-B-79) 


where  C°  = dH°/dT  is  the  molar  constant  pressure  heat  capacity  of  the 

tt,  pl 

i gaseous  constituent  in  its  standard  state, 


l <r 

i P 


■ C 


Z cn 

'pi  j=p+apij  pi*(3> 


and 


, / _ 


<*  Mfir  £p 

0 dT  p 


(I  -B-80) 


(I-B-81) 


The  quantity  Re^  is  the  contribution  to  the  constant  volume  heat  capacity 
of  a mixture  of  weight  M0  resulting  from  gas  phase  imperfections. 

Summing  the  Gibbs  free  energy  of  the  condensed  and  gaseous  phases 
gives  the  equation  for  the  specific  Gibbs  free  energy  of  the  mixture  as 


Mog 


£ s ji*  + ,r 
J=i  3 3 1=1 


Vi 


Q-B-82) 


An  alternative  expression  for  g in  terms  of  the  stoichiometric  constant 
q is  obtained  as 


Mog 


+ J=P+iVi*(j) 


s 

+ Z n ; u 
t=i  it  i 


P 


* 


c 

j=p+l^ijV(j) 


(I-B-83) 


by  multiplying  Eq . (I-B-57)  by  p.*  and  Eq.  (I-B-58)  by  ^ and 

adding  the  resulting  expression  to  Eq.  (I-B-82).  Combining  Gibbs’  con- 
dition for  chemical  equilibrium  Eq.  (I-B-47a)  with  Eq.  (I-B-83)  gives 
the  free  energy  of  the  equilibrium  mixture  as 


Mog  = 


Z q p + 
)=i  J J 


c 

J=P+i  J i*<j) 


( I -B  -84) 


I -B  -30 


since  the  bracketed  term  in  Eq.  (I-B-83)  vanishes.  Combining  Eq.  (I-B-84) 
with  Eqs.  (I-B-28)  , (I-B-SO),  and  (I-B-S3^  gives  the  equation  for 
g('"1  , ....  - ) in  terms  of  *-  as 

1 C+2 


Mog/TRT  = 


c 

T q r + 
j=i  J 3 


h q - + i 

j=p+i  j c+l  J=p+l 


nr 


(I-B-85) 


The  specific  internal  energy  of  the  mixture  is  calculated  with 
the  thermodynamic  identity 


e(’-  , ....  n ) = h - pv  (I -B-86) 

a c+a 

and  specific  entropy  is  calculated  with  the  identity 

s(  - J = (h  - g)/T  (I-B-87) 

l c+a 

c.  State  Conditions 

The  state  conditions  expressed  in  Eq.  (I-B-69)  by  F - F = 0 

c+i  c+a 

define  which  one  of  the  seven  pairs  of  state  variables  <Po.To)»  (Po^ho)  , 
CpOtSo).  (Po  .v0)  , (vc  ,T0)  , (v0  ,e0)  , or  <v0  ,s0>  is  specified  at  the 
thermodynamic  point  of  interest.  The  case  when  the  pressure  pq  is 
specified  will  be  considered  before  the  case  when  the  specific  volume 
v0  is  specified. 

When  the  pressure  of  the  system  is  specified  as  pq  , is 

defined  by  the  equation 

F = in  pr  - in  p = 0 (l  -B-88) 

c+x 

and  p is  calculated  with  the  gas  phase  equation  of  state.  At  (pq ,T0) 

points,  the  condition  + = constant  is  used  to  specific  that  T = T0 

c+s 


I-B-31 


5>  :.4S  9SS&$%& 


and  consequently  the  equation  1' 


0 is  not  required.  The  iteration 


procedure  used  to  calculate  the  equilibrium  composition  is  simplified 
because  the  rank  of  the  Newton-fiaphson  matrix  is  reduced  to  c + 1 and 


the  condition  r 


c+a  ‘c+2 


=0  is  imposed  at  each  iteration.  At  (po ,bc) 


and  (PqiS01  points,  Fc+-  is  defined,  respectively,  by  the  equations. 


F = M°<HC  “ h>/RT  = 0 


(I-B-89) 


F , * Mo(s0  - (h  - g)/T)R  = 0 


(I-B-90) 


where  h is  evaluated  with  Eq.  (I-B-78)  and  p = Po  and  g is  evaluated 
with  Eq.  (I-B-85) . At  (Po.Vo)  points  F is  defined  by  the  equation 

q+2. 


F = PoMoCvq  - v)/RT  = 0 


(I-B-91) 


and  v is  evaluated  with  Eq.  (I-B-75) . 

When  the  specific  volume  vc  is  specified  and  the  pressure  is 
not,  then  F^^  is  defined  by  the  equation 


F - pMo  ( v0  - v)/RT  = 0 

c+i 


( I -B  -92) 


and  v is  evaluated  with  Eq.  (I-B-75).  At  a (v0 ,T0)  point,  the  tem- 
perature T0  is  specified  in  the  same  way  as  for  the  (po  ,T0)  point. 

At  (v0 ,e0)  and  (v0,s0)  points,  F is  defined  respectively  by  the 

C + 2 

equations 


F = Mo(e0  - h + pv0)/BT  * 0 

C + 2 


(I  -B-93) 


F = Mq(s0  - (h  - g)A)R  = 0 

C + 2 


( I -B -94) 


where  h is  evaluated  with  Eq,  (I-B-78)  and  g is  evaluated  with 
Eq  . (I  -B -85)  . 


I-B-32 


d.  Partial  Derivatives  with  Respect  to  th<'  Iteration  Parameters 


It  is  convenient  for  notational  purposes  to  expand  the  expressions 
(I-B-65)  in  a linear  Taylor  series  about  an  approximation  (r^ , n°+a) 

to  the  solution 


(rt  , ) and  write 

i c+a 


c+a  _ 

,Er  ,4  = FO 

k*i  jk  k j 


J = 1 , . . , , c + 2 


(I-B-95) 


to 


with  = - dF^/dr^,  = ^ “ 11°,  and  the  superscript  0 used 

denote  that  the  quantity  so  designated  is  evaluated  at  the  approximation 
to  the  state  cf  the  system.  The  equations  for  the  partial  derivatives 
of  the  state  variables,  obtained  in  the  expressions  for  D by  differ- 


entiating  F partially  with  respect  to 


according  to  the  chain  rule, 


will  be  presented  before  the  expressions  for  D 

jk 

(1)  Partial  Derivatives  of  the  Mole  Numbers 

Equations  for  the  partial  derivatives  of  the  mole  numbers 
with  respect  to  the  iteration  parameters  are  derived  with  Eq.  (I-B-29) 
and  (I-B-55).  Differentiating  Eqs.  (I-B-29)  and  (I-B-55)  with  respect  to 
r gives  the  equations 


and 


d*. 


i dr. 


(I-B-96) 
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- + T.  a o 

r~i  dn  Sr  J+1  pij 
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A 

with  k - 1,  c + 2,  the  partial  derivatives  dp/dr^  and  ST/d^ 


Jk 


(I-B-97) 
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ft 

given  by  the 

equations 
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and  a defined  by  the  equation 

pi,c+a 


i,c+a  B" 


c+a , 


dg. 


dT 


(I-B-99) 


since  is  a function  of  temperature  only.  It  follows  from  Eq.  (I-B-33) 

that  1$^  can  be  expressed  in  terms  of  the  reduced  enthalpies  of  the 

gaseous  species  in  their  standard  states  as 


^i,c+a  *i  j=p+i  ^i;j*i*(j) 


(I -B-100) 


Setting  dn  /Br,  = n bf  /dr-  and  > = ST  /Bn  in  Eq.  (I-B-97)  and 

r It  r^r  k lr  i r 

rearranging  terms  leads  to  the  following  equations  with  i = 1,  ....  s: 
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(I-B-103) 


where  the  derivatives  of  the  activity  coefficient  are  given  according 

to  Eq.  (I-B-31)  by  the  equations. 


ar  P .a  * 

i _ !*  Mq  Bp  dp 

Bn^  J0  RT p Bn^Bn^  p 
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(I-B-104) 


(I-B-105) 


and 


I-B-34 


di  2 A 
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(I-B-106) 


Although  the  integral  expression  for  dT\/d  in  T is  given  here,  this 
derivative  can  be  readily  obtained  bv  differentiating  the  explicit 
expression  for  I*  directly. 

The  relationships  in  Eqs.  (I-B-101)  to  (14-103)  constitute 

c + 2 sets  of  s simultaneous  linear  equations  for  and  can  be 

solved  for  these  derivatives  when  the  mole  numbers  associated  with  a set 

(rj,  , r>  ) of  the  iteration  parameters  have  been  calculated  with 

* C+2 

Eq.  (I-B-71)  in  the  inner  iteration.  The  corresponding  derivatives  of  the 
gaseous  mole  numbers  can  then  be  calculated  with  Eq.  ( I -B-96) . 


The  Eqs.  (I-B-101)  to  (I-B-103)  become  simplified  for  an 
ideal  gas  mixture  with  =0,  i =1,  ...,  s.  In  this  case,  the  matrix 
of  Eqs.  (I-B-101)  to  (I-B-103)  reduces  to  the  unit  matrix,  and  the 


equations  have  the  solution 


i * 1 , * . . • s 

k — 1 , « . • , c + 2 


(I-B-107) 


The  abbreviation 


n dr- 
i k 


i = 1, 
k - 1, 


, . , c +2 


(I-B-108) 


is  introduced  for  convenience  in  writing  other  derivatives  '-f  the  state 
with  respect  to  the  iteration  parameters.  Differentiation  of  Eq.  (I-B-57) 
gives  the  partial  derivatives  of  the  mole  numbers  of  the  condensed  species 


- r a C n 
i“i"ij  ik  i 


3 1 » » • » , Pi 

k = 1 c + 2 


I-B-35 


and  differentiation  of  Eq.  (l-B-18)  with  p - p'  gives  the  derivatives 
of  the  total  number  of  gaseous  Moles  n as 


_ c n 

i*i**i,c+l  ik  i 


k * 1,  ....  c + 2 (I-B-110) 


with  « . defined  by  Eq.  (I  HB-54)  . 

''i , c+l 


(2)  Partial  Derivatives  of  the  Pressure 

It  is  convenient  to  derive  equations  for  these  partial 
derivatives  with  Eq,  (I-B-74)  for  the  pressure  rewritten  as 


in  p * An  ♦ + An  n + n 

‘c+l 

Differentiation  of  Eq.  (I-B-lll)  with 
following  equations 


8 An  p 

8 to  ♦ 1 

Bn 

dTV 

Bn,  n 

k 

8 An  p 

B An  * ( 1 

Bn 

&n  . 

‘c+l 

Bn  n 

C+i 

Bn 

‘c+i 

8 An  p 

8 fa  ♦ 1 

Bn 

Bn 

C+2 

Bn  n 

C+2 

Bn 

C+2 

(I-B-lll) 

respect  to  then  leads  to  the 

k = 1,  ....  c (I-B-112) 

(I-B-113) 

(I-B-114) 


where  4 is  considered  to  be  a function  of  tl  , . . . , r . Differentiating 

t c+a 

♦ by  the  chain  rule  and  asking  use  of  Eqs.  (I-B-98)  and  (l-B-108)  leads 
to  the  following  expressions  for  the  partial  derivatives  of  ♦ 
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(I  -B-1I6) 
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d in  $ _ 5 ia  $ d In  $ + $ In  $ 
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c+s  i 


Equations  (I-B-112)  to  (I-B-114)  for  the  partial  derivatives  of  the  pres- 
sure are  evaluated  with  the  aid  of  Eq.  (I-B-llO)  for  the  derivatives  of 
n and  Eqs.  <I-B-115)  to  (I-B-117)  for  the  derivatives  of  $. 

(3)  Partial  Derivatives  for  the  Condensed  Phases 

The  thermodynamic  properties  of  a condensed  phase  are  assumed 
to  be  explicit  functions  of  temperature  and  pressure  specified  by  its 
complete  equation  of  state.  They  can  therefore  be  regarded  as  implicit 
functions  of  the  iteration  parameters  because  of  the  implicit  dependence 
of  the  pressure  on  these  parameters.  The  partial  derivatives  of  the 
thermodynamic  properties  of  a condensed  phase  with  respect  to  n are 
derived  because  they  are  required  to  calculate  the  partial  derivatives 
of  the  mixture . 


The  partial  derivatives  of  the  molar  volume  V*  and  the 

* th  J 

molar  enthalpy  H of  the  J condensed  component  are  expressed  in 
J 

terms  of  the  identities  contained  in  Eqs.  CIHB-36) , (I-B-37),  end  (I-B-40) 
with  j = 1 , . . . , p as 

a d An  p 
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and 
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k = 1,  ....  c + 1 (I-B-120) 
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I-B-37 


The  partial  derivatives  of  the  reduced  chemical  potential  $t*/ST  of  the 

th  J 

j condensed  component  arc  oxpressed  in  terms  of  the  identities  in  Eqs. 

(I-U-40)  and  (l-U-41)  with  j • 1,  ....  c as 
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P.  <P.T> 

, 

RT 


* d in  p 

tA 


k=l,  2,  ...,c+l  (I -B-122) 
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RT 


J df) 


c+s 


(US-123) 


In  a point  calculation  with  the  pressure  specified  to  be 
Po • the  properties  of  the  condensed  phases  can  be  evaluated  in  each  itera- 
tion with  po  instead  of  with  the  current  value  of  the  pressure  calculated 
with  the  gaseous  equation  of  state.  In  this  case,  the  properties  of  the 
condensed  phases  can  be  regarded  as  a function  of  temperature  only,  and 

the  equations  for  the  partial  derivatives  of  V*.  H*.  and  U*/BT  are 

J J 

obtained  by  setting  d An  p/  3tl  * 0 <k  = 1 c + 2)  in  Eqs.  (I-B-118) 

to  (I-B-123)  . 


(4)  Partial  Derivatives  of  the  Specific  Volume  of  the  Mixture 


The  partial  derivatives  of  the  specific  volume  with  respect 
to  the  iteration  parameters  can  be  written  as 
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(l-B-125) 


(I-B-126) 


I -B -38 


These  equations  can  be  rewritten  with  Eqs.  (I-B-118)  and  (I-B-119)  as 
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K - 1,  2,  ...  , c (I-B-127) 
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The  Eqs.  (I-B-127)  to  (I-B-I29)  are  used  to  evaluate  the 

derivatives  of  F and  F at  a point  where  the  specific  volume  of 
c+i  c+a 

the  mixture  is  specified  to  be  v0 . 

(5)  Partial  Derivatives  of  the  Specific  Enthalpy  of  the  Mixture 
The  identity 


1 

o An  T 


(I-B-130) 


together  with  Eqs.  (I-B-120)  and  (I-B-121)  for  the  derivatives  of  H*  , 

J 

Eq,  (I-B-110)  for  the  derivatives  of  n , and  Eq.  (I-B-79)  for  the  frozen 
heat  capacity,  are  used  to  express  the  partial  derivatives  of  the  specific 


enthalpy  as 
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for  k = 1,  2,  c and 
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The  equations  for  the  derivatives  3h/&r^  at  a point  where  the  pressure  is 

specified  to  be  po  can  be  obtained  by  setting  & An  p/&t^  =0  (k  = 1 , 

....  c + 2)  in  Eqs.  (I-B-131)  to  (I-B-132)  because  in  this  case  H*  is 

3 

a function  of  tenperature  only. 

(6)  Partial  Derivatives  of  the  Specific  Free  Energy  of  the 
Mixture 

The  derivatives  d (Mo are  readily  obtained  by 
differentiating  Eq.  (I-B-85)  as 


r^“  * q.  for  k = l,  2 c 
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(I-B-135) 


(I-B-136) 


I-B-40 


(7)  Partial  Derivatives  of  the  Specific  Energy  and  Specific 


Entropy  of  the  Mixture 

The  equations  for  the  derivatives  and  fts/Cir 

are  obtained  by  differentiating  the  identity  for  e Eq.  (I-B-86)  and 
the  identity  for  s Eq.  (I-B-87)  and  substituting  the  equations  for  the 
resulting  derivatives  that  have  already  been  presented. 

(8)  Equations  for  the  D Derivetives 


The  partial  derivatives  of  the  state  variables  witn  respect 

to  the  iteration  parameters  presented  in  the  previous  paragraphs  can  now 

be  used  to  formulate  the  equations  for  the  D , derivatives.  The  equa- 

jk 

tions  for  the  condensed  components  (j  = 1,  . . . , p)  will  be  presented 
first,  then  the  equations  obtained  from  the  stoichiometric  conditions 
<J  = p + 1,  . . . , c) , and  then  those  obtained  from  the  state  conditions 
(j  = c + 1 , c + 2)  . 


The  equations  for  D with  j = l , ...,  p are  obtained 

Jk 


a *** 

D =6 - i 

3 k Jk  Br  RT 
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k - 1 , > * . , c + 2 


( I -B-137) 


by  differentiating  Eq.  <I-B-67)  partially  with  respect  to  r-  . Combining 
Eq.  (I-B-137)  with  Eqs.  (I-B-122)  and  ^I-B-123)  gives  the  equations 
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which  reduce  to  the  equations 


Jk 


k = 1,  2 c + 1 (l-B-140) 
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J ,C  + 2 


(X-B-141) 


for  the  point  calculations  when  the  pressure  of  the  system  is  specified 
to  be  p.;  . 


obtained  as 


The  equations  for  B with  j = p + 1,  c are 

J*4 


dn. 

_s  i 


Djk  = A $ij  k = P + 1>  ••**  c + 2 (I-B-142) 

by  differentiating  Eq.  (1-8-68) , and  they  are  rewritten  as 

V = ii  SlJClk"l  k » P + 1 c + 2 (I-B-143) 


with  the  Eq.  (I-B-108). 

The  equations  used  in  the  seven  different  point  calculations 
for  °jk  With  J“c+1»c+2  are  obtained  by  differentiating  the 
state  conditions.  The  case  when  the  pressure  of  the  system  is  specified 
to  be  pq  will  again  be  considered  before  the  cose  when  the  specific 
volume  is  specified  to  be  v0 . 


obtained  as 


At  points  where  p = p0 , the  equations  for  D are 

c+i  ,k 


D 


c+i  ,k 


d in  p 

dr. 


k = l,  2,  . , , c + 2 (I-B-144) 


by  differentiating  Eq.  (I-B-88)  partially  with  respect  to  r , and  the 

k 

derivatives  are  evaluated  with  Eqs.  (I-B-112)  to  (I-B-114),  At  (p0  ,T0> 
points,  the  D derivatives  are  not  required  because  the  condition 

C+3  ,k 

^C+E  - 0 is  employed.  The  equations  for  D ^ ^ are  obtained  at  (p0  ,ho) 
points  as 
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D 


c+a  ,k 


Mg.  '3h_ 


RT  o 


k = 1 , 2,  ....  c + 1 (I-B-145) 


fc 


and 


c+a ,c+a 


C-r-a 


(I-B-146) 


by  differentiating  Eq.  (I-B-89);  at  <Po(So)  points  as 

k = 1 , 2 , . ..,  c + 1 (I-B-147) 


It  = “s.  'dh  ) 5 Mpg 

C+2  ,k  RT  &t^  RT 


and 


C+2 , C+2 


= “2.  f{*!L  ) . h1  . __L_  «2l 

RT  L\3r  /o  J 5n  RT 


(I-B-148) 


c+a 


c+a 


by  differentiating  Eq.  (I-B-90) ; and  at  (Po ,v0)  points  as 


d = BbMsl/Sl.} 

c+a,k  RT  Vdn^o 


k-1,2,  . ...  c 1 (I-B-149> 


and 


D 


c+a ,c+a 


- as.  n-|!_)  t (V0  . 

RT  lA  d+  /0 
C+8 


V) 


(I-B-1^0) 


by  differentiating  Eq.  (1-45-91)  . The  (dh/3’-,  )0  and  Qv/3r'.  )o  deriva- 

k k 


tives  are  evaluated  with  the  equations  obtained  by  setting  O £n  p/3ru ) = 0 

k 


in  Eqs.  (I-B-131)  to  (I-B-133)  and  in  Eqs . (I-B-127)  to  (I-B-129)  , and 
the  ddiog/RTJ/Sr^  derivatives  are  evaluated  with  Eqs.  (I-B-134)  to 


(I -B-136) . 


At  points  where  v * vq  and  the  pressure  is  not  specified, 


the  equations  for  D are  obtained  as 

c+i  , k 


„ pMo  r^v  , x 3 Xn  p'l 

De« ,k  - mi*-- ^ - v)  "sr“J 


A 


k = 1 c + 1 (I-B-151) 
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aaaaaa.  -- — — 1 1 


and 

D 

c+i ,c+a 


p**c  rav 

rt  La- 

c+a 


+ <v0 


a in  p\ 
6r  / 


c+a 


(I-B-152) 


by  differentiating  Eq,  (I-B-92) , and  the  derivatives  are  evaluated  with 

Eqs.  (I-B-112)  to  (I -B -114)  and  Eqs.  (I-B-127)  to  (I-B-129)  . At  {v0  ,T0) 

points,  the  D . derivatives  are  not  required  because  the  condition 
c+a,k 

D ~ 0 is  employed.  The  D are  obtained  at  (vn  .eo)  points  as 

C+2  e+a  ,k 


D 

c+a  ,k 


Mo  T3h 

rt  la-^ 


pv0 


5 in  p~! 

a-  J 

It 


k = 1,  . , c + 1 (I-B-153) 


and 


Tj 

c+a , c+a 


■c+a 


S p , 

pv0  + (e0 


Sn 


C+2 


'1 

h + pvo)^ 


(US-154) 


by  differentiating  Eq.  (I-B-93) , 


D _ Mq  Sh  5 Mpg 

c+a,k  RT  dr,  S-  RT 
k k 


and  at 


(v0 ,s0)  points  as 
k = 1 , c + 1 


and 


D 

c+a ,c+a 


Mo  /d h 
RT  \br 

C+s 


- h 


d Mpg\ 

9r  RT  J 
c+a 


(1  -B-155) 


(I-B-156) 


by  differentiating  Eq.  (I-B-94) , 

e'  The  Equilibrium  Partial  Derivatives 

Equations  will  be  derived  for  a complete  set  of  equilibrium 
first-order  partial  derivatives  of  the  system  so  that  the  equations  for 


i 

! 


I 

f 


■( 

| 

> 

i 


1 


i 

\ 
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- — — - * — ngr 


any  other  equilibrium  partial  derivative  can  be  obtained  with  thermodynamic 
identities. 


A convenient  set  of  complete  derivatives  for  use  in  hydrodynamic 
calculations  consists  of  the  equilibrium  specific  heat  capacity  at  constant 
volume 

'de\  ,'Bh\  fbp\ 

iSOV 

V V V 


«w  - (I)  -(g)-  (I) 


(I-B-157) 


and  the  following  derivatives  of  the  equilibrium  (e-p-v)  equation  of 
state 

Q 

T /As\  d /& T\ 

“(— ) -1  (I-B-158) 


O = 


VBe\  = T /ds\  _ x & _p  / £T\ 

p\Bv/  p \Bv/  p VBv/ 


- _ i _ T /as\ 

B ~ v \Bp/  ~ v \3p/ 


[v  /BTN 

kBp/  v \3p/  v V Bp/ 

v v \ 


(I-B-159) 


The  adiabatic  exponent  x used  hydrodynamic  calculations  is  related 
to  a and  ^ by  the  identities 


/BInp\ 
“ U in  W 


a + i 


J8 


(I-B-160) 


It  is  convenient  to  rewrite  the  Eqs.  (I-B-157)  to  (I-B-159)  and  change  the 


independent  variables  to  r and  ri  with  the  identities 

C+J,  C+a 


(R  ■ (**-)  Aw-) 

T 'c+l  - • '«**i  . 


( I-B-161) 


Tl 


C+l 


C-f  s 


c+a 


1 (I?).,  ■ & - kb)  & <— > 


’c+a 
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c+a 


c+l 


c+i 


c+a 


c+i 


c+a 
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so  that  the  quantities  C^,  a,  and  £ can  be  evaluated  with  the  partial 
derivatives  presented  in  Section  d.  Care  must  be  taken  however  to  re mem- 
ber that  (d/5-  ) and  (5/dr 


(d/5-  ) 

c+X  r 


) 

c+a  r 


ore  equilibrium  partial  deriva- 


C+  S?  ' C+  )> 

tives  related  to  the  frozen  derivatives  5/&r^  (k  = 1 , . . . , c + 2) 

by  the  equations 


0/5-  ) 

c+i  , 


= J-  (d/dr.  ) (ar  /dn  ) 
k=x  k k c+i  . 


a/5r 


c+i 


(I-B-163) 
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(a/a-  > 

c+a 


I O/ari  H3r  /an  ) 

k=i  k k c+a  . 


+ a/ar 


c+a 


(I-B-164) 


c+ 1 


C+i 


and  that  the  equilibrium  derivatives  of  the  iteration  parameters  +^(k  = 1» 


....  c)  are  required  for  their  evaluation.  The  derivatives  On. /5+  ) 

k c+x  r 


denote  the  rates  of  change  of  the  parameters  with  respect  to 


c+a 


c+i 


when 


c+a 


is  kept  constant  and  the  system  remains  in  chemical 


equilibrium.  Similarly,  the  On.  /dr>  ) 

V ‘c+a  r, 


derivatives  denote  the 


rates  of  change  of  the  +v  parameters  C+1  with  respect  to  r> 

k c+a 


when  kept  constant  and  the  system  remains  in  chemical  equilibrium. 


The  equations  for  C , a,  and  g are  rewritten  as 

V ** 


MqC 


v _ Mo  Oh\  /McPV\  r_/g  in  pt  1 

r ■ r vs)  ■ vtr)  J 
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and  the  identities  {I-B-160  and  (I-B-162)  are  used  to  express  the  partial 

derivatives  appearing  in  these  equations  in  terms  of  r and  n as 

'c+t  c+s 


RT 

MoP 


5 P\  /^MqP  Sv  \ 

dri  / / \ RT  Sr  / 

c+i  n c+i  n 

c+a 
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(US-168) 
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c+a 


The  application  of  Eqs.  (I-B-163)  and  (I-B-164)  to  the  pressure,  the  specific 
volume,  and  the  enthalpy  then  gives  the  equations  relating  the  equilibrium 
partial  derivatives  to  the  frozen  derivatives  presented  in  Section  d. 

The  equations  for  the  derivatives  of  the  pressure  are  written  as 


/5  In  p\ 

\ Sr  / 


c+i  r 


c+a 


P H + —V, 

Sr  k=l 

C+i 
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and 
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„c  5 In  p 
k=l  Srk 
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(I-B-172) 


(I-B-173) 
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with  the  frozen  derivatives  of  the  pressure  given  by  the  Eqs.  (I -B-112) 
to  tI-B-114).  The  equations  for  the  derivatives  of  the  specific  volume 
of  the  mixture  and  written  as 


MOP  j 
RT  1 

f dv  \ 

1 = “a£  9v  , rc  Sv  i 

f 3rlk\ 
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(I-B-17S) 


with  the  frozen  derivatives  given  by  the  Eqs.  (I-B-127)  to  (I-B-129), 
and  those  for  the  specific  enthalpy  as 


“2-1 

( 

dh  + ^c  Mj, 

RT  ’ 

^ dr  J RT 

dr  k=l  RT 

Vd-  ) 

C+l  r 

c+a 

'c+i 
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and 
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(I  -B-17  7) 


with  the  frozen  derivatives  given  by  Eqs.  (I-B-131)  to  (I-B-133) . The 
equations  for  the  equilibrium  derivatives  of  the  iteration  parameters 


’Vtk  = 1,  ....  c),  (B+  ^ ) , and  (dr,  /dr.  ) 


k c+i 


c+s 


k c+a'  c+i 


are  obtained 


for  j = 1,  c by  differentiating  Eqs.  (I-B-67)  and  (I-B-68)  and 

setting  dFj/d’Tj,.  = - DJk  as 


kil 


= D 


'c+a 


J .c+i 


(1 -B-17  8) 


and 
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£*  D . a) 
k=i  jk  Tt  c+a 


= D 


C+i 


J.c+S 


(I-B-179) 


expressions  for  D for  the  condensed  components  j = 1,  p are 

JK 

given  by  Eqs.  (I-B-138)  anl  (i-B-139)  , and  expressions  for  D for  the 

)* 

gaseous  components  3 = p + 1,  ....  c are  given  by  £q.  (I-B-143) . 
Equations  (I-B-178)  and  (I-B~179)  are  solved  for  the  equilibrium  deriva- 
tives with  the  values  of  D obtained  by  solving  Eq.  (I-B-95)  for  the 

jK 

equilibrium  composition  with  the  Newton-Raphson  method. 


Equations  (I-B-165)  to  (I-B-179)  are  used  to  evaluate  the  complete 

> at  a point  where  the  system  has  attained 
chemical  equilibrium  without  a condensed  component  changing  phase. 


set  of  derivatives  < Cy>  a, 


f • Phase  Change  of  a Condensed  Component 

The  foregoing  treatment  of  condensed  species  cannot  be  used  to 

solve  the  equilibrium  problem  when  the  condensed  species  are  considered 

capable  of  existing  in  different  phases.  The  treatment  must  be  extended 

to  determine  when  a condensed  species  is  present  as  a single  phase  or  as 

an  equilibrium  mixture  of  the  phases.  The  important  practical  case  of 

melting  when  the  phases  are  solid  and  liquid  is  presented  to  illustrate 

the  basic  method  used  to  solve  this  type  of  problem. 

Gibbs'  conditions  for  equilibrium  used  previously  to  test  the 

a priori  assumption  made  about  the  presence  of  condensed  constituents 

must  now  be  used  to  determine  the  relative  stabilities  of  the  liquid  and 

solid  phases  of  a condensed  component.  Let  u*  and  u*  denote  the 

chemical  potentials  of  the  solid  and  liquid  phases  of  the  jth  condensed 

component.  N*  and  N*  enote  their  respective  mole  numbers,  T*  « 
sj  ij  J 

T*(p)  denote  the  Clsusius-Clapeyron  equations  for  the  variation  of 
J 

melting  points  with  pressure,  and  f^  denote  the  fraction  a* 


I-B-49 


liquid  present  at  the  melting  point.  Gibbs'  conditions  for  equilibrium 
of  the  solid  and  liquid  phases  can  then  be  written  for  j = 1 P as 


p*  (p.T  < T*>  < p*  <p,T  < T*) 
•J  J *3  3 

^(P.T  > T.)  < ^<P,T  > T*) 

* 4‘P-T  ■ 


(I-B-180) 

(I-B-181) 

(l-B-182) 


with  J =1#  . ...  p.  Equation  (I-B-180)  gives  the  stability  condition  for 

the  solid,  Eq.  (I-B-181)  the  stability  condition  for  the  liquid,  and 

Eq.  (I-B-182)  the  stability  condition  for  a mixture  of  solid  and  liquid 

at  a melting  point.  In  an  equilibrium  mixture  the  3th  condensed 

component  will  exist  as  a solid  with  = 0 when  Eq.  ( I -8-180)  is 

satisfied,  as  a liquid  with  f = 1 when  Eq.  (I-B-181)  is  satisfied, 

and  as  a mixture  of  solid  and  liquid  with  0 < f £ 1 when  Eq.  (I-B-182) 

A 

is  satisfied. 


The  Eqs.  (I-B-180)  to  (I -B-182)  apply  to  all  condensed  components, 
but  at  this  stage  melting  of  only  one  of  them,  the  3th  , will  be  considered. 
It  is  then  convenient  to  write  the  following  equations  for  the  properties 
of  the  j1*1  condensed  component  constituent. 


Of  = f fit.  + (l  - f.)af 
3 Jt  Jtj  i sj 
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(I  -B-188) 
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Equations  ( I -B -1 83)  to  (I-B-189)  are  used  with  those  already  presented  for 
the  aixture  to  calculate  the  equilibrium  thermodynamic  state  at  a specified 
point  when  the  jtto  condensed  component  can  occur  as  a liquid,  a solid,  or 
a mixture  of  liquid  and  solid  at  the  melting  point.  As  in  the  previous 
treatment  of  the  equilibrium!  problem  with  condensed  species, the  validity 
of  the  calculated  state  is  determined  by  its  compatibility  with  the 
assumptions  required  to  perform  the  calculations.  It  should  be  remem- 
bered, however,  that  the  question  whether  a condensed  species  is  present 
as  a liquid  or  a solid  is  not  a problem  in  a (p,T)  calculation  because 
the  phase  is  determined  by  the  specified  values  of  p and  T and  the 
Claus ius-Clapeyron  equation.  In  the  particular  case  when  the  free  ener- 
gies of  the  solid  and  liquid  are  found  to  be  equal,  the  (p,T)  calculation 

is  performed  at  the  solid  phase  boundary  with  f - 0. 

*/ 

Consider  the  other  six  point  calculations  when  the  phase  assump- 
tion is  fixed  by  the  values  of  T and  p,  say  (T0  .pt,)  f chosen 
to  start  the  iterative  procedure  to  determine  the  thermodynamic  state. 

Only  the  case  when  the  choice  of  the  solid  as  a component  is  valid  and 

* * 

N ^0  will  be  considered.  Suppose  further  that  ip  (T0 iPo)  > 0 so 
sj  J 

that  the  calculation  is  perfortaed  with  f „ - 0,  and  let  T and  p 

l c c 

denote  the  values  of  temperature  and  pressure  resulting  from  the  calcu- 
lation. Since  the  point  calculation  was  performed  with  f - 0,  it  is 

Jbr 

necessary  to  check  that  the  calculated  state  is  compatible  with  this 
phase  assumption.  The  calculated  state  is  an  equilibrium  state  when 


I -B-51 


either  y <T  ,p  ) > 0 end  the  solil  is  acre  stable  then  the  liquid,  or 
^ j c c 

to i (T  ,p  ) * 0 end  the  solid  is  is  equilibrium  with  the  liquid  at  the 
J c c 

♦ 

solid  phase  boundary.  But  it  is  not  an  equilibriua  state  when  au  (T  ,p  ) < 0 

3 c c 

and  the  liquid  is  more  stable  than  the  solid,  for  then  the  phase  assumption 

and  the  calculated  state  are  incompatible.  When  au*(T  ,p  ) <0,  values 

j c c 

of  To  and  p-  are  chosen  so  that  &y*(T0,po)  < 0 and  the  point  calcu- 

* 

lation  is  performed  again  with  f = 1.  Kow  with  a 0,  the  calculated 

X XJ 

state  ia  an  equilibrium  state  when  Ay*<T  ,p  ) < 0 and  the  liquid  is 

3 C C 

more  stable  than  the  solid  or  is  in  equilibriua  with  the  solid  at  the 

liquid  phase  boundary,  but  not  when  Ap*(T  ,p  ) > 0 and  the  solid  is 

J c c 

more  stable  than  the  liquid.  When  au*(T  ,p  ) >0,  it  follows  that 

“j  c c 

Ap*  =o  at  the  point  of  interest , and  the  thermodynamic  state  must  be 
J 

in  the  mixed  phase  region  with  f in  the  range  0 < f < 1. 

Z if 


When  *0,  the  value  of  t . in  the  mixed  phase  region  must 

be  determined  to  perform  a point  calculation  except  at  a (p.T)  point 

where  T*  = T*(p)  and  the  value  of  f can  be  chosen  at  will.  Since 

3 3 * 

Ay  < 0 in  the  first  calculation  with  f = 0,  and  Ap*  >0  in  the 
3 * J 

second  calculation  with  f ^ = 1 , it  is  convenient  to  use  the  equilibrium 
* 

condition  Ay  = 0 to  determine  the  equilibriua  value  of  f . Since  f 
3 X A 

lies  in  the  range  0 < f^  < l,  a bounded  linear  approximation  is  used  in 

the  iteration  scheme  to  determine  the  value  of  f that  satisfies  the 

X 

equation  Ay*  = 0. 


The  order  of  the  procedure  must  obviously  be  changed  when  the 
3th  condensed  component  is  first  assumed  to  be  present  as  a liquid  and 
not  as  a solid.  The  development  of  a general  procedure  for  solving  the 
equilibrium  problem  when  more  than  one  of  the  condensed  species  can  melt 
will  be  based  on  the  treatment  for  one  species  presented  in  this  section. 
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g.  Determination  of  a Mixed  Phase  Li  tie 

The  routine  for  perforating  mixed  phase  line  calculations  is 

* * 

called  MIXXED.  The  mixed  phase  line  T = T (p)  of  a condensed  species 

J <3 

in  the  (p-T)  plane  is  determined  by  solving  the  following  equation 

V*  (T,p)  = P*  (T,p)  (I-B-190) 

sj  Xj 

expressing  the  equilibrium  condition  for  the  solid  and  liquid  phases,  if 
values  of  temperature  are  specified,  then  Eq.  (I-B-190)  is  solved  with  a 
Newton-Raphson  iteration  scheme  to  determine  the  corresponding  values  of 
pressure.  Similarly,  if  values  of  pressure  are  specified,  Eq.  (I-B-190) 
is  solved  to  determine  the  corresponding  values  of  temperature. 


The  thermodynamic  state  in  a mixture  containing  a melting  com- 
ponent can  be  calculated  readily  when  the  mixed  phase  line  of  this  compo- 
nent has  been  determined.  The  calculation  is  performed  at  a (p,T)  point 

on  the  mixed  phase  line  with  a specified  value  of  f . However,  calcu- 

X 

lations  of  only  two  of  these  equilibrium  states,  namely,  those  with 

f =0  and  f = 1 , have  been  programmed  into  the  code.  These  values 
X X 

of  f were  chosen  purposely  to  define  the  phase  boundaries  and  the 
mixed  phase  region  of  the  thermodynamic  system  in  the  pressure-volume 
(p-v)  plane. 


h.  glmse  Ctwggg  of  Condensed  Species 


The  treatment  of  melting  developed  for  one  species  was  extended 
to  t - s species  to  obtain  a more  realistic  treatment  of  thermodynamic 
systems  containing  multiple  condensed  phases.  The  equilibrium  conditions 
expressed  by  Sqs.  (I-B-180)  to  (I-B-182)  must  now  be  used  to  test  the 
phase  assumptions  for  t - s species  rather  than  for  one.  A calculated 
state  is  an  equilibrium  state  when  it  satisfies  the  phase  assumptions 
required  to  make  the  calculation.  As  discussed  in  I -B-4f , however,  there 
is  no  problem  in  (p,T)  calculations  because  the  phases  are  determined 
by  the  specified  values  of  p and  T and  by  the  Clausius -Clapeyron 
equations  for  the  species. 

Consider  the  other  six  point  calculations  when  the  phase  assump- 
tions are  fixed  by  the  values  of  T0  and  Pq  chosen  to  start  the  itera- 
tive procedure  to  determine  the  thermodynamic  state.  It  is  convenient  to 
label  those  species  set  as  solids  with  m = 1 , . . . , a and  those  set  as 
liquids  with  m = a + 1 , ...,  a'  5 t - s,  so  that  the  point  calculation 

is  performed  with  f =0  for  m * l „ ...,  a and  with  i.  -I  for 

Xm  Lm 

m = a + 1,  It  follows  from  Section  I -B-4f  that  a calculated 

a 

state  is  an  equilibrium  state  when  it  satisfies  the  conditions  (T  ,p  ) a 

u m c c 

0 for  m * l , . . . , a , and  &u  (T  ,p  ) SO  for  m = a + 1 , ...,a/. 

^m  c c 

When  any  one  of  these  conditions  is  not  satisfied,  the  calculated  state 
is  not  an  equilibrium  state  and  it  is  necessary  to  choose  new  values  of 
T0  and  pq , say  (To'.Po*),  and  repeat  the  point  calculation  with  the 
new  phase  assumptions  determined  by  the  Clausius -Clapeyron  equations. 

The  values  of  T0 ' and  ' are  based  on  the  results  of  the  last  calcu- 
lation and  on  the  fact  that  the  phase  line  of  each  species  with  an  incor- 
rect phase  assumption  intersects  the  line  segment  joining  (p0 ,T0)  and 

(p  ,T  ) in  the  (p  - T)  plane.  These  points  of  intersection  are  first 
c c 
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determined , and  the  values  of  Tc ' and  Pg * are  chosen  to  change  the 

phase  assumption  of  the  particular  species  whose  intersection  point  lies 

closest  to  <ps.T0).  Thus  the  phase  assumption  of  only  one  species  is 

changed  by  this  procedure,  and  (po  ,T0)  lies  on  one  side  of  its  phase 

line  and  (p^ # ,T0 lies  on  the  other.  Calculation  of  the  state  now 

gives  new  values  of  T and  p , say  (T*  ,p  ')  , and  the  procedure 

c c c c 

is  continued  until  the  calculated  state  agrees  with  the  phase  assumptions 

used  to  perform  the  calculation.  If  the  phase  assumptions  for  a particu- 

lar  species,  say  with  m » b,  are  found  to  be  incorrect  in  two  successive 

point  calculations,  the  state  must  lie  in  the  mixed  phase  region  of  this 

species  with  0 < f 1.  In  this  case,  the  phases  of  the  other  con- 
J6D 

densed  species  are  fixed,  and  the  thermodynamic  state  in  the  mixed  phase 
region  is  calculated  using  the  method  presented  previously  in  Section 
I -B-4f . 

i , The  Equilibrium  Partial  Derivatives  in  a Mixed  Phase  Region 
This  section  will  be  completed  at  a later  date. 
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I-C.  THEORETICAL  BASIS  FOR  NONEQtS  ILIBRIUM  CALCULATIONS 
IN  PARTIALLY  FROZEN  SYSTEMS 


1 . Introduction 

The  theoretical  basis  for  making  nonequilibrium  calculations  is  an 
extension  of  that  developed  by  R . S,  Brinkley,  Jr.  for  making  equilibrium 
calculations,*  Methods  are  presented  here  for  calculating  the  thermodynamic 
state  of  a nonideal  heterogeneous  system  in  partial  equilibrium  when  the 
pair  of  state  variables  can  be  chosen  from  the  following  set:  [(p,T) , 

(p,h) , <p,s) , (v,T)  , (v,e)  , (p,v) , and  (v,s)].  The  nonequilibrium  states 
considered  are  those  in  mechanical  and  thermal  equilibrium,  with  part  of 
the  composition  frozen  and  part  of  it  in  chemical  equilibrium.  Frozen 
species  are  defined  as  species  with  mole  numbers  that  are  fixed  and  do 
not  change  in  calculations  of  the  thermodynamic  state.  But  frozen  con- 
densed species  are  allowed  to  change  phase  according  to  their  equations 
of  state  as  in  equilibrium  calculations.  In  calculating  the  thermodynamic 
state  of  such  a system,  the  mole  numbers  of  the  frozen  species  must  be 
specified,  and  the  mole  numbers  of  the  other  species  must  be  calculated 
using  the  equilibrium  conditions  and  the  stoichiometric  conditions  given 
in  Section  I-B  of  the  TIGER  documentation.  It  is  important  to  note 
however  that  the  choice  of  the  frozen  composition  may  not  be  compatible 
with  the  total  stoichiometry  of  the  system.  The  frozen  species  should 
be  chosen  so  that  the  mass  balance  equation  for  the  system  can  be  satis- 
fied. In  other  words,  it  is  necessary  to  ensure  that  there  is  a feasible 
solution  to  the  mass  balance  equations  for  the  nonf rozen  species.  This 
problem  was  treated  in  a way  similar  to  that  used  for  the  linear  programming 
* 

Presented  earlier  in  Section  I-B-4. 

I-C-l 


problem,  end  a routine  called  FEASIB  was  written  to  determine  the  choice 
of  frozen  constituents  compatible  with  the  total  stoichiometry  of  the 
system. 


It  is  important  to  remember  that  the  description  of  the  composition 
presented  in  Sections  I-B-2  and  I-B-3a  applies  equally  well  to  nonequi- 
librlum  system  es  to  equilibrium  system  because  it  was  formulated  with 
the  conservation  of  mass,  which  is  valid  for  ail  thermodynamic  systems. 

In  fact,  it  was  for  this  reason  that  the  notation  for  labelling  frozen 
constituents  was  introduced  in  Section  I -B  rather  than  in  I-C,  The  dis- 
cussion of  equations  of  state  presented  in  I-B-3b  is  also  valid  for  our 
system  in  partial  equilibrium  because  the  pressure  and  temperature  are 
well  defined  by  the  assumptions  of  mechanical  and  thermal  equilibrium. 

For  convenience  in  presentation,  the  treatment  given  in  this  section 
of  systems  in  partial  equilibrium  will  parallel  that  given  in  section 
I -B  for  systems  in  chemical  equilibrium.  Consequently,  equations  that 
are  the  same  as  those  in  Section  I-B  will  be  repeated  in  special  cases 
only. 


2,  Calculations  with  a Partially  Frozen  Composition 

■ Equilibrium  Conditions  and  Iteration  Parameters 

When  the  system  is  assumed  to  be  in  partial  equilibrium,  the 
frozen  composition  is  chosen  and  the  remainder  is  calculated  with  the 
equilibrium  and  stoichiometric  conditions.  These  conditions  will  be 
expressed  in  terms  of  the  notation  for  frozen  constituents  introduced 
in  Section  I-B-3a.  The  condensed  constituents  assumed  to  be  present  and 
in  equilibrium  with  the  gaseous  phase  are  labelled  with  j = 1,  2,  ...,  p. 
These  constituents  will  always  be  chosen  as  components.  The  condensed  con- 
stituents assumed  to  be  frozen  and  not  in  chemical  equilibrium  are  labelled 
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with  m = p •*■  1 , . . . , p“,  and  the  remainder  that  are  assumed  to  be 
absent  (but  to  satisfy  the  equilibrium  conditions)  with  p'+  1,  ...,  t - s. 
The  gaseous  species  assumed  to  be  components  and  in  equilibrium  with  the 
rest  of  the  system  are  labelled  with  g = 1,  ...,c-p,  the  remaining 
equilibrium  species  are  labelled  with  g * c - p + 1,  . s*,  and  those 
assumed  to  be  frozen  with  g=s,+l,  ...,  s. 

The  stoichiometric  conditions  for  the  mixture  are  written  with 
the  i/(m)  subscript  introduced  on  page  I-B-12,  and  with  terms  summed 
from  p + 1 to  p/  to  account  for  frozen  condensed  species  as 


“ « ,,  v N + £ fl  n + N = q 

m=p+lpi/(m)j  m i=r  ij  i j j 

/ 

P s 

ZB,,  . N •*■  £ fl  n = q 

m=p+l*i  (ra)j  m i=rij  i 


3 = 1,  2 


j = p + 1,  . . . , c 


a-c-r. 


(I-C-2) 


The  conditions  for  the  chemical  potentials  of  gaseous  constituents  are 
written  as 


Hf(g) 


P 

E. 

J=l 


^f<g> 


c 

T 


j-p+i^{(g)  A*(j> 


g = i, 


(I-C-3) 


and 

^i(g)  ^ j=l^i(g) 


j=p^l^(g>J^(j) 


K — s 1 + 1,  . . . , s 


(I-C-4) 


The  first  of  these  expressions  is  used  to  obtain  equilibrium  conditions, 
and  the  second  is  included  to  account  for  the  fact  that  s - s'  gaseous 
species  are  frozen  and  are  not  in  chemical  equilibrium  with  the  remainder 
of  the  system.  Remembering  that  Eq.  (I-C-3)  contains  c-p  identities 
for  the  constituents  chosen  as  components,  the  equations  for  the  chemical 
potentials  of  the  species  in  our  partially  frozen  mixture  can  be  written 
as : 


^f(g) 


=l5f(g) 


* 

P,  + 


J J 


j=p+l^f(g) j^fX j) 


B = 1 » 


(I  -C-5a) 


I -C-3 


,y. -y-y'y  - 


Pi(g)  ^ J«l^t(g) j^j  +j=p+A(g)jtlf(j) 


g = s'  + 1,  . ..,  s (I-€-5b) 


1<J) 


J * 1, 


Cl-C-Sc) 


V«  = M '<»>/.  + d4+A  WfO>  ' £ • = P + >’  •••• p/  <^'5d) 

»VU>  = M'(m>/0  + < < B " p'  + 1 1 * 3 (lHC‘5e) 

The  gaseous  constituents  in  equilibrium  must  satisfy  (I-C-5a) , and  those 
that  are  not  in  equilibrium  must  satisfy  (I-C~5b).  The  condensed  con- 
stituents assumed  a priori  to  be  present  and  in  equilibrium  must  satisfy 
(I-C-ic),  the  frozen  condensed  species  must  satisfy  (I-C-5d),  and  those 
assumed  to  be  abseil  must  satisfy  (I-C-5e)  . Equations  (I-C-l)  , (l-C-2)  , 
and  (I-C-S)  are  sufficient  for  calculating  the  equilibrium  composition 
in  the  mixture  at  one  of  our  specified  state  points  when  the  mole  numbers 
of  the  frozen  species  (N^  £0,  ra=p+l,  , ...  p'  and  n^  , g = s ' + 1 , 

. ..,  s)  have  been  properly  assigned  and  a complete  equation  of  state  of 
the  mixture  is  known. 


Since  u.  is  known  as  a function  of  the  n.  's , the  s - c + p 

independent  equilibrium  conditions  together  with  the  c stoichiometric 

conditions  constitute  a set  of  s'  + p equations  for  calculating  the 

s#  + p unknown  equilibrium  mole  numbers.  With  regard  to  the  previous 

treatment,  note  that  the  number  of  equilibrium  conditions  and  the  number 

of  unknown  mole  numbers  for  the  gaseous  species  are  both  reduced  by 

s - s'  by  the  inclusion  of  s - s'  frozen  species.  Here  again,  the 

computed  equilibrium  composition  must  be  checked  to  test  the  validity 

of  the  a priori  assumptions  made  about  the  presence  (N  SO,  j = 1 , ....  p) 

J 

and  the  absence  (N  = 0,  m = p + 1,  . . . , t - s)  of  condensed  constituents, 
m 
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The  assumptions  are  valid  when  the  computed  equilibrium  composition  is 


y 

\ 


found  to  satisfy  the  condition  N 20  and  that  expressed  bv  Eq.  (I-C-5e). 

J 

Most  of  the  equations  for  calculating  the  equilibrium  gaseous 
composition  (Section  I-B-4a)  can  be  used  for  our  system  in  partial  equi- 
librium when  the  constituent  indices  are  changed  to  account  for  the 
presence  of  frozen  species.  It  is  necessary  to  change  the  index  i with 
values  i = 1,  ....  s to  1(g),  with  g = 1,  . ..,  s',  and  to  change 
the  values  of  m to  range  from  m = p + 1,  . t - s to  m = p'  + 1, 

....  t - s.  There  is  no  need  to  change  the  j index,  however,  because 
it  refers  to  components. 


Thus  the  equations  for  . = in  n.,  . (g  = 1 s')  are 

i(g)  i(g) 

obtained  by  substituting  i(g)  for  i in  Eqs.  (I-B-48),  (I-B-49), 

(I-B-55) , (I-B-60)  and  (I-B-61) . The  expressions  for  g*  are  obtained 

m 

by  substituting  p'  for  p in  Eqs.  (I-B-51),  (I-B-52)  , and  (I-B-56)  . 
The  equations  (I-B-50),  (I-B-53),  (I-B-54),  and  (I-B-59)  used  to  define 


n. 


^i,c+l’ 


and 


remain  the  same.  Similarly,  the  equations  for 


'U  g) 


and  j are  obtained  by  changing  i to  f(g)  in  Eqs , (I-B-62) 


and  (I -B-63) , but  Eq.  (I-B-64)  for  is  unchanged. 


The  equation  modifications  required  to  make  the  iterative  pro- 
cedures for  calculating  the  equilibrium  composition  applicable  to  a 
system  in  partial  equilibrium  will  now  be  discussed.  Recall  that  Brinkley's 
approach  to  the  equilibrium  problem  was  to  construct  a set  of  c + 2 
equations 


F ( 

J 


l • 


C + 2 


) = 0 


(I-C-6) 


so  that  the  set  of  *•  satisfying  these  equation  defines  the  equilibrium 
composition.  The  definitions  of  for  the  system  in  partial  equilibrium 

are  as  follows: 
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x*  /Rt  - n =o 


J - 1 (I-C-7) 


’ = q - E J}  n + - fl/  N 

J J ij  i m=p+l  i(ra)j  m 


J = P + 1,  ....  c (I-C-8) 


F =0 
c+i 


F =0 
c+a 


(i-e-9) 


These  equations  differ  from  those  used  for  the  system  in  chemical  equi- 
librium only  by  the  summation  terra  for  the  frozen  condensed  species  in 

Eq.  (I-C-8)  . The  iterative  equations  for  F are  the  same  as  those 

J 

presented  in  (I-B-66) , but  the  following  equations  for  their  derivatives, 
corresponding  to  those  given  in  (I-B-70) , 


8F  s 

^ * - Iff 


*{(*> 


g=lPf  (g)  j "f  (g) 


J * P + 1 c (I-C-10) 


must  be  modified  with  the  conditions, 


= s'  + 1,  ....  s (I-C-ll) 


to  ensure  that  the  concentrations  of  the  frozen  gaseous  species  do  not 
change  when  the  composition  of  the  part  of  the  system  in  equilibrium  is 
calculated.  The  equations  for  ’ corresponding  to  those  presented 

for  in  (I-B-71)  to  (I-B-73) , are  as  follows: 


Wi(g)^i(l)’  •'•’**<■'))  ' " *l(g)  " r£(g) 


+ ->(^1 J = 0 

l(g)  C+» 

g — 1,  **«,  s 


(I-C-12) 


' dw  » x 
r=l  iU> 


’ ^i(s))^f(r)  " W{(g)(^i<l)’  *‘”  ^f(s')) 


g - 1,  ....  s' 


(I-C  ’ 1) 
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^i(r) 


^(g) 
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, /^(g)N 

C{(g)f(r)  Vdr;t(i);  t(r> 


(I-C-14) 


' l^5( 


3n^  .SOi.  .< 
i(g)  t(r) 


b.  State  Equations 

The  state  equations  are  used  to  calculate  the  thermodynamic 
state  of  our  system  in  partial  equilibrium  when  its  equilibrium  composi- 
tion has  been  determined.  Only  the  equations  differing  from  those 
presented  in  Section  I-B-4b  will  be  given. 

The  equation  for  the  specific  volume  of  the  mixture  is  written 


) = 4 + £ N V*  + £ S V*) 

2 P M-o  \ J=1  j 3 m=p+l  m m / 


a-c-15) 


and  the  corresponding  equation  for  the  specific  enthalpy  at 

P * P#  s n A 

Moh  - Z ,N  H + I N H*  + Z n + plio/O  + RT(e  - n)  (l-t-16) 

J=1  3 3 m=p+l  m m i=l  i i ' 


with  the  imperfection  term  defined  by  Eq.  (I-B-77).  The  alternative 

expression  for  h(r^ , ) is 

C + 2 


P * £ * £ 

Moh  - 77  N H*  + Z N H + Z ,1 

j=l  j j m=p+l  m m j=p+l 


/ S 

Hi*(j)^qj  m=p+ A '("O 
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+ RT  Z,0.  n.  + pMo/p  + RT( f - n) 

i=l  ljC+s  i 


and  the  equation  for  the  frozen  heat  capacity  of  the  mixture  is 
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J=1  3 PJ  m=p+l  m pm  j=p+l  pi*(j) 


> \ j m=p+l^i(m)  j'W 


(I-C-18) 


Z <A.C°)n  + R(e'  - n) 

1=1  i p 1 T 


Equations  (1-8-82)  through  (I-B-85)  for  the  specific  free  energy  must  be 
altered  because  part  of  the  composition  is  frozen.  Summing  the  Gibbs' 
free  energy  of  the  condensed  and  gaseous  phases  of  the  mixture  gives  the 
specific  free  energy  of  the  mixture  as 
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(I-C-19) 


Following  a procedure  with  Eq.  (I-C-19)  similar  to  that  presented  for 
Eq.  (1-8-82)  gives  the  following  expressions  for  the  free  energy 


Mog  = 


3=llVi  3 + J=p+l^i(  j)'*j 


? N fp*  - 

ia=p+l  m\  m 


q-  + — " jMi'(m) /j  " j=p+l8i'(»)/ni) 


g=s/+lRf(g)(tAt(g)  j t^(g) 


(I-C-20) 
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3 " J=P+l8i\g)/i*(j). 


Mog/RT  = 


£ q r;  + Z a (r  + p“  % /RT) 
J=1  3 3 3=P+1  3 c+1 
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g=I'+inf(g)^f(g)/R'  j‘^l^f(g)  3*3  J=pfl^{(g)  j^i+Ij)^^ 


The  last  two  summation  terms  in  Eqs.  (I-C-20)  and  (I-C-21)  are  for  the 
frozen  species,  and  they  disappear  if  the  system  is  allowed  to  come  to 
equilibrium. 
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c*  State  Conditions 

The  state  conditions  expressed  in  (I-C-9)  by  F = F =0 

c+1  c+2 

define  which  one  of  the  seven  pairs  of  state  variables  (po >To) , (Poj^c) 

(Po.So)*  (Po  ,va)  , (v0,T0),  (va  ,eo>  • or  Cv0 ,s©)  is  specified  in  the 

thermodynamic  calculation.  The  equations  used  to  define  F and 

c+1 

F c+2  for  these  seven  different  point  calculations  are  obviously  the  same 
as  those  already  presented  in  Section  I-B~4c  of  the  TIGER  documentation 
and  will  therefore  not  be  given  here. 

d.  Partial  Derivatives  with  Respect  to  the  Iteration  Parameters 

Partial  derivatives  with  respect  to  the  iteration  parameters 
are  required  to  evaluate  the  derivatives  F = - D ^ in  the  equation 

\ = F°  (I-C-22) 

k j 


which  are  used  to  solve  the  set  F ^ = 0(j  = 1,  . . . , c + 2)  with  the 
Newton-Raphson  technique  in  the  TIGER  code. 


(1)  Partial  Derivatives  of  the  Mole  Numbers 

The  equations  for  these  partial  derivatives  are  derived 
from  the  relationships 


*f<g)  = ln  Dt(g) 


g = 1,  ....  s'  ( I -C  -23) 


*!(g)  ~ " gf(g) 


c*l. 
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g = 1 s'  (I  -C-24) 


by  the  procedure  given  in  Section  I-B-4d-l  of  the  TIGER  documentation. 

They  are  therefore  essentially  the  same  as  the  equations  given  in  I-B-4d-l 

and  can  be  readily  obtained  from  them  by  changing  i to  i(g)  , with 

/ 

g = 1.  s',  and  r to  t(r)  , with  r = 1,  ...,  s',  where  necessary. 
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fa 
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and  including  the  additional  conditions  for  the  frozen  species,  ft  is 
necessary  to  change  the  i index  in  Eqs . (l-B-96),  (l-B-97),  (1-0-99) 
through  (I-B-103),  (I-B-107),  and  (l-H-IOK),  the  r index  in  Kqs . 
(I-B-101)  through  Eq.  (1-6-103) , and  to  add  the  following  equations  for 
frozen  species, 
dn. 


f(g) 

r* 

k 


= o 


C*tg>,k  “ ° 


g = s'  + 1,  ...,s  (I-C-25) 

g = s'  + l , . . . , s 


k - 1 , . . . , c + 2 


( I -C  -26) 


as 

m 


= 0 


ra  - p + 1 , . . . , p 

, , n Cl-C-27) 

k = 1 , . . . , c + 2 


It  should  be  noted  that  although  the  subscript  i(g)  could  be  used  in 
(I-B-110) , there  is  no  need  to  change  the  notation  because  the  summation 
is  taken  over  all  the  gaseous  species. 


(2)  Partial  Derivatives  of  the  Pressure 

The  equations  (I-B-lll)  through  (I-B-117)  for  the  pressure 
derivatives  derived  previously  for  systems  in  chemical  equilibrium  are 
also  used  for  systems  in  partial  equilibrium. 


(3)  Partial  Derivatives  for  the  Condensed  Phases 


The  partial  derivatives  of  the  molar  volume,  molar  enthalpy, 

and  chemical  potential  of  the  condensed  species  in  a partially  frozen 

system  are  obtained  by  adding  the  derivatives  of  the  frozen  species  to 

those  given  for  the  equilibrium  species  in  (l-B-118)  through  (I-B-123). 

The  equations  for  the  frozen  species  are  generated  iron  Eq.  (l-B-118) 

through  Eq.  (i-B-123)  by  changing  j(j  =1,  . . . , p)  to  m(m  = p + 1, 

...  p') . The  equations  for  V*  are  given  as  an  example. 
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V*{p,T) 

RT  -V  m 
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«*«  %i=-S 

m n 

k 


k 1,2,  ....  c + \ (t_c_28) 

m “ p + 1 , . . . , p 


* -f-  v*(,.d  . e«f  |-^ 

RT  or  o no  a a d-r 


in  = p + 1 p'  (r-C-29) 


<4)  Partial  Derivatives  of  the  Specific  Voltaae  of  the  Mixture 


The  derivatives  of  the  specific  volume  are  obtained  by 
differentiating  Eq.  (I-C-15)  with  respect  to  the  iteration  parameters. 
They  differ  from  the  derivatives  for  the  equilibrium  mixture  by  a sun- 
Buitiuu  term  over  the  condensed  frozen  species  as  shown,  for  example,  in 
the  following  equations  corresponding  to  Eqs.  (I-B-127)  through  (I-B-129) 


es& -J- . Vs-e ( S»  rv  * * ,» sv) 

RT  j=l  6r\^  v.  J*1  jhj  j m=p+l  ® m m/ 


(I-C-30) 
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-*  ( £,N  jSV  + — .N 

„ \ J=1  j J J m=p+l  m in  W 
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J+1  J J J m=p+l  m » m RTp 


(I  -C  -32) 


(5)  Partial  Derivatives  of  the  Specific  Enthalpy  of  the  Mixture 

Here  again,  the  partial  derivatives  of  the  specific  enthalpy 
of  the  partially  frozen  mixture  differ  from  those  for  the  equilibrium  mix- 
ture by  summation  terms  for  the  frozen  species.  The  equations  are  readily 
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tfcvmi  «***^*i* WS» * 


obtained  by  adding  - 9 in  to  (I-B-131),  by 

adding  - d in  p/d*  .[  ? _K  OrtS'*')  to  Eq.  {I-B-132)  and  by  adding 
c+I\B— p+1  m m m/ 

- 9 in  pArs  ( 1 M oftflf)  to  Eq.  (I-B-133). 

c+2nb=p+1  st  m w 

(6)  Partial  Derivatives  of  the  Specific  Free  Energy  of  the 
Mixture 

The  derivatives  5(Mog/RT)/9r^  are  obtained  by  differen- 
tiating (I-C-21)  as 

* Mqk  jg'  /9  in  p \ . 5 /^jjg>  n \ 

9r^k  RT  qk  m=p+l  b^  9t^  m itm)k/  g=s'+l  i(g>\  9^  ik/ 

k = 1 c 


= ! q + ? j(. 

J=P+1  J B=p+1  Bt\ 


d in  p 


^ ^Km)e+1 


s Cg)  \ 

+ g=s '+inf(g)Van  , _ Cg) c+1/ 
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)*"  RT  i=p+l  iXi*( j)  np=p+l  ha  9r  j j=p+l  Km)  j*i*( j) / 


v *(g)  o v o' 

+ g=s/+lnt(g)\9'-  *t(g)  + j=p+l*f(g)jXi*(j)< 

c+z 


e.  The  Nonequi librium  Partial  Derivatives 


The  identities  for  the  complete  set  of  partial  derivatives 
(C  , a,  #>  given  in  Section  I-B-4e  of  the  TK&R  documentation  apply 
to  any  system  whose  thermodynamic  state  can  be  described  by  two 


I-C-12 


independent  variables.  Consequently  the  equations  given  in  I -8-4e  in 
terns  of  the  iteration  parameters  can  be  used  to  evaluate  these  deriva- 
tives for  the  system  in  partial  equilibrium  when  care  is  taken  to  account 
for  the  fact  that  part  of  the  composition  is  frozen.  This  fact  is 
automatically  taken  account  of  in  TIGER , however,  because  the  derivatives 

used  in  the  evaluation  are  them- 


'W,  , ““ 

*c+l 


c+2 

selves  evaluated  with  the  values  of  C..  calculated  for  the  equilioriura 


ik 


species  and  with  the  values  of  C = 0 for  the  frozen  species. 


f . Phase  Changes 


The  treatment  of  phase  changes  for  systems  in  partial  equili- 
brium is  exactly  the  same  as  that  presented  for  equilibrium  systems  in 
Sections  I -B -4f  through  I -B-4h  of  the  TIGER  documentation. 
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I-D.  THEORETICAL  BASIS  FOR  NONEQU I LIBRIUM  CALCULATIONS 
IN  COMPLETELY  FROZEN  SYSTEMS 

1 . Introduction 

The  method  for  computing  the  thermodynamic  state  of  a nonideal 
heterogeneous  system  with  a frozen  composition  will  be  considered  as  a 
limiting  case  of  the  method  described  in  Section  I-C  for  computing  the 
state  in  systems  with  a partially  fi.zen  composition.  Here  again  the 
nonequilibrium  states  are  considered  to  be  in  mechanical  and  thermal 
equilibrium  and  the  condensed  species  are  allowed  to  change  phase  as  in 
equilibrium  calculations,  but  the  mole  numbers  of  all  the  species  are 
fixed,  in  any  one  of  the  seven  point  calculations  for  determining  the 
thermodynamic  state,  the  composition  is  chosen  to  satisfy  the  stoichio- 
metric conditions,  and  the  calculation  is  performed  with  the  concentrations 
of  the  species  fixed. 

2.  Calculations  with  a Completely  Frozen  Composition 


The  completely  frozen  system  is  described  in  terms  of  the  indices 
introduced  to  label  species  in  Section  I-B-3  of  the  TIGER  documentation 
by  the  conditions  p = s7  = 0 and  p7  = t - s.  Consequently,  the  equa- 
tions used  to  perform  point  calculations  for  completely  frozen  systems 
are  readily  obtained  by  setting  p=s7  = 0 and  p7  = t - s in  the 
equations  presented  in  Section  I-C-2  for  the  partially  frozen  system. 

For  this  reason  they  will  not  be  given  here.  The  equations  for  C , Ot, 

v 

and  0 are  the  same  as  those  presented  in  Section  I-B-4e  of  the  docu- 
mentation for  equilibrium  systems.  In  this  case,  however,  the  deriva- 
tives are  evaluated  with  C =0  for  i - 1,  ....  s and  k - 1 , .... 

ik 

c + 2 because  all  the  species  are  frozen. 

I -D-l 
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I -E . INCORPORATION  OF  AN  EQUATION  OF  STATE 
INTO  TIGER 


1 , Gaseous  Equation  of  State 

The  incorporation  of  a nonideal  gaseous  equation  of  state  into  TIGER 
to  perform  the  calculations  discussed  in  the  previous  sections  is  by  no 
means  a trivial  task.  It  is  necessary  to  use  the  imperfection  term  4> 
to  derive  expressions  for  the  thermodynamic  quantities  required  to  per- 
form the  calculations.  These  quantities  include  the  frozen  partial 
derivatives  of  $,  b in  4t/b  in  T,  b in  <$/d  in  p,  and  d in  ^/Sn^  ; the 
imperfection  integrals  ( and  ; and  the  activity  coefficients  T.  and 
their  frozen  partial  derivatives.  Equations  for  these  quantities  are 
derived  with  the  Becker,  Kistiakowsky , Wilson  (BKW)  relationship  to 
demonstrate  the  incorporation  of  a nonideal  equation  of  state  into  the 
code.  The  corresponding  expressions  for  the  JCZ2  and  JCZ3  equations  of 
state  developed  by  Jacobs,  Cowperthwaite , and  Zwisler  are  well  documented 
in  the  final  report,  "improvement  and  Modification  to  TIGER  Code"  written 
under  Contract  N 69921 -72-C -0013,  and  will  therefore  not  be  presen*  t here. 

The  form  of  the  Becker  equation  of  state  proposed  by  Halford, 
Kistiakowsky,  and  Wilson1  and  modified  by  Fickett  and  Cowan  has  become 
known  a rj  the  BKW  equation  of  state.  It  is  written  here  in  terms  of  the 
universal  constants  OL,  and  Q as 

p = nRTpfft'Mo  . 

<3>  - <ftx)  = 1 + xe^  (I-E-l) 
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X = k£/Mc(T  + 6)a 
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(I-E-2) 


win-re  * is  a scale  factor  Introduced  for  computat tonal  convenience, 

* t h 

kk^  denotes  the  covolurae  cons  tart  of  the  i gaseous  constituent,  and 
k denotes  the  covolurae  of  the  mixture.  Differentiating  Eq.  (I-E-l) 
leads  to  the  following  equations  for  the  frozen  derivatives  of  $ , 
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(I-E-6) 


Let 


g 


denote  the  internal  energy  of  the  gas  mixture  contained 


in  a unit  mass  of  mixture,  and  let  denote  the  molar  internal  energy 

of  the  i**1  constituent  in  its  standard  state  taken  as  the  hypothetical 
ideal  gas  at  unit  pressure.  Then  the  imperfection  integral  ( for  the 
gas  mixtur^  can  be  written  as 
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nr  ^ f.n  n*l  nA 

(I-E-7) 
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and  the  integration  performed  along  an  isotherm  with  the  composition  fixed 

♦ 

as  that  of  the  mixture. 


Performing  the  integration  in  Eq.  (I-E-7)  with  Eqs.  (I-E-l)  to  (I-E-3) 
gives  the  imperfection  integral  £ for  the  BKW  equation  of  state  as 
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nog . 

T + 8 


(I-E-9) 


and  the  corresponding  imperfection  integral  is  obtained  by  differ- 

entiating (Eq.  I-E-9)  partially  with  respect  to  T as 
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The  expression  for  the  activity  coefficient  is  derived  with  the  following 
equation  introduced  earlier 
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Performing  the  integration  in  Eq . (l-E-il)  with  Eqs,  (I-E-l)  , (I-E-2), 
and  (I-E-5)  gives  the  activity  coefficient  I\  for  the  BKW  equation  of 
state  as 
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( I -E  -12) 


And  differentiation  of  Eq . (i-E-12)  gives  the  corresponding  equations  for 
the  partial  derivatives  of  needed  for  tne  equilibrium  calculation  as 

* See  Bef . 3 for  a careful  discussion  of  nonideal  mixtures  and  the 
standard  state. 
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The  thermodynamic  identities 


i _ Mqp  a in  p 
a in  4 RTp  dn 


&C  x _ Mpp  f*  6 in  p~ 

a An  P RT$  L & in  T- 


a in  T 


are  useful  fcr  checking  the  expressions  derived  for  T and  ( with  a 

i 

A 

particular  p = p(p,  T,  nx , . ..,  n ) equation  of  state. 

s 

2 . Condensed  Equations  of  State 

Since  routines  for  condensed  equations  of  state  are  still  being 
developed,  this  section  will  be  completed  at  a later  date. 
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I HF  , ONE -DIMENSIONAL  DETONATION 
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l . C ha proan-J ouguet  Theory  of  Detonation 
a.  Brief  History 

The  one-dimensional  model  of  a detonation  wave  was  formulated 
by  D.  Chapman,1  and  independently  by  E,  Jouguetf  at  the  turn  of  the  century. 
Its  purpose  was  to  account  for  wave  propagation  in  combustible  gases  at 
speeds  of  the  order  of  2000  meters  per  second.  Since  experiments  had 
shown  the  propagation  to  be  supersonic,  the  front  of  the  detonation  wave 
was  assumed  to  be  a reactive  shock  discontinuity.  The  exothermic  process 
in  the  discontinuity  was  viewed  as  maintaining  the  velocity  of  the  wave, 
and  a thermodynamic-hydrodynamic  theory  based  on  this  assumption  was 
developed  to  calculate  that  velocity.  The  theory  combines  the  energy 
yield  of  the  combustion  reaction  with  the  conservation  laws  of  mass, 
momentum  and  energy.  These  well  established  laws,  together  with  the 
equation  of  stnte  and  the  Chapman-J ouguet  hypothesis,  define  the  problem 
and  allow  completion  of  the  calculations. 

Calculated  velocities  in  gaseous  explosive^  are  in  good  agree- 
ment with  experiment,  which  is  partly  to  be  explained  by  the  fact  that 
the  perfect  equation  of  state  can  be  validly  applied.  Successful 

4 

calculations  of  detonation  velocity  in  gases  led  Schmidt  to  apply  the 
Chapman-J ouguet  theory  to  condensed  explosives.  Here,  however,  an 
equation  of  state  at  several  hundred  kilobars  is  required,  but  in  general 
vio  satisfactory  equations  are  available  in  this  pressure  range. 


I-F-l 


b.  Rankine-Rugoniot  Relationships 


The  Rankine-Hugoniet  relationships  express  the  continuity  of  | 

mass,  momentum,  and  energy  across  a shock  discontinuity.  Since  the  rela-  f 

S I 

t ions hips  are  well  known,  they  will  be  presented  without  derivation.  Let  i 

£ 

D,  u,  and  p = 1/v  denote  shock  velocity,  particle  velocity,  and  f 

i 

density,  and  let  the  subscripts  0 and  1 denote  values  at  the  bottom  - 

and  top  of  tho  shock  discontinuity.  Then  the  Rankine-Hugoniot  relation-  / 

ships  for  a shock  pi'opagating  at  velocity  D in  a stationary  material 


(t*D  = 0) 

can 

be 

written  as 

AoD 
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Pi(D  - ux) 

(I-F-l) 

Pi 

- Po 

= 

(I-F-2) 

Pi«i 
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jOfeDtej  ~ + \ u*) 

(I-F-3) 

Equation 

(I  -F 

-1) 

expresses  the  conservation  of  mass,  Eq. 

( I -F  -2)  expresses 

the  balance  of  momentum,  and  Eq.  (I-F-3)  expresses  the  first  law  of 
thermodynamics  under  the  assumption  that  the  shock  process  is  adiabatic. 
Moreover,  since  the  shock  process  is  irreversible,  the  second  law  of 
thermodynamics  is  expressed  by  the  following  inequality  for  the  specific 
entropy,  Si(ei,Pi)  > s0 (e0  ,pq)  . Combination  of  Eqs.  (I-F-l)  and  (I-F-2) 
gives  the  equation 

px  - Pc  » (D/v0)s(v0  - vx)  <I-F-4> 

for  a straight  line  in  the  (p  - v)  plane  that  passes  through  the  initial 
conditions  (Pc,vo)  and  is  called  the  Rayleigh  line.  Combination  of 
Eqs.  (I-F-l),  (I-F-2)  and  (l-F-3)  gives  the  equation 

ej  - e0  = | (Pi  + Po)  (v0  - Vi)  (I-F-5) 


% 


I-F-2 


which  relates  the  thermodynamic  variables  across  the  shock  and  is  called 
the  Hugoniot  equation.  An  initial  condition  <e3  ,po  ,v0)  , an  (e-p-v) 
equation  of  state  of  shocked  material,  and  the  Hugoniot  equation  define 
a curve  in  the  (p-v)  plane  called  the  Hugoniot  curve  centered  at  (po,v0). 
The  compressive  portion  of  the  Hugoniot  curve  specifies  the  locus  of  states 
attainable  from  an  initial  state  <ee ,Po »v0)  by  shocks  with  different 
velocities.  For  a nonreactive  shock,  the  Hugoniot  curve  passes  through 
its  center  point  <Po >vo> , but  this  is  not  so  for  a reactive  shock 
because  of  the  energy  change  in  the  reaction  process.  In  calculating 
Hugoniot  curves  with  Eq.  (I-f-5),  the  heat  of  reaction  is  automatically 
accounted  for  by  the  difference  in  the  specific  internal  energies  of  the 
standard  states  of  the  reactants  and  the  products  of  combustion. 

For  exothermic  waves,  the  compressive  part  of  the  Hugoniot 
curve  is  called  the  detonation  branch  of  the  Hugoniot,  and  the  expansive 
part  is  called  the  deflagration  branch.  Since  a shocked  state  must  satis- 
fy Eqs . (I-F-4)  and  (I-F-5),  the  intersection  of  the  Hugoniot  curve  cen- 
tered at  (po,v0)  and  the  Rayleigh  line  of  slope  -(D/v0)3  passing 
through  (po ,v0)  defines  the  thermodynamic  state  (ej ,pj ,vi)  behind 
either  a nonreactive  or  resetive  shock  discontinuity  traveling  with  con 
stant  velocity  D in  stationary  material  with  pressure  Po  and  specific 
volume  v0 . 

Since  the  slope  of  the  Rayleigh  line  must  be  positive  to  satisfy 
the  conservation  of  mass,  and  momentum  hiwever,  the  detonation  branch  of 
the  Hugoniot  curve  is  terminated  at  the  point  where  ex  = e©  and  = 

v0  , and  the  deflagration  branch  of  the  Hugoniot  curve  is  terminated  at 
the  point  where  = ho  and  Pi  “ Po  • These  termination  points  are  for 
obvious  reasons  called  the  constant  volume  explosion  point,  and  the  con- 
stant pressure  combustion  point,  and  will  be  denoted  by  the  superscripts 
e and  c.  The  pressure  at  the  constant  volume  explosion  point  can  then 
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be  written  formally  as  p = p<e0 ,v0) , and  the  volume  at  the  constant 

e 

pressure  combustion  point  as  v = vth^po). 

c . Chagan^jougaet  Hypothesis 

Since  the  states  on  a Hugoniot  curve  centered  at  (pq ,vc) 
represent  the  locus  of  states  connected  to  an  initial  state  (eo ,po ,v0) 
by  shocks  traveling  at  different  velocities,  an  additional  condition  is 
needed  to  determine  a unique  detonation  velocity.  In  other  words,  the 
three  conservation  laws,  together  with  a complete  equation  of  state  for 
detonation  products  in  chemical  equilibrium,  are  insufficient  to  calcu- 
late the  variables  D,  %,  px  , fa  , ej  , which  characterize  the  detonation 
process.  Since  there  are  five  unknowns,  but  only  four  equations  to 
calculate  them,  an  additional  equation  is  needed  to  complete  the  theory. 
The  fifth  equation,  known  as  the  Chapman-douguet  condition 

D a + cx  (I  -F  -6} 

where  cx  is  the  velocity  of  sound  in  the  detonation  products  at  the  top 
of  the  shock  discontinuity,  follows  directly  from  the  Chapman-douguet 
hypothesis  that  the  head  of  a rarefaction  wave  at  the  shock  discontinuity 
travels  at  the  same  speed  as  the  shock  itself. 

It  follows  from  the  Chapman-douguet  condition  that  the  stable 
detonation  wave  is  represented  in  the  (p-v)  plane  by  the  Rayleigh  line 
(see  Figure  1-1)  through  (po  ,v0)  that  has  a point  of  tangency  with  the 
Hugoniot  curve.  This  point  of  tangency  defines  the  compressed  state  at 
the  front  of  the  detonation  wave  and  is  called  the  Chapman-douguet  point. 
It  is  of  interest  to  derive  these  properties  of  the  Chapman-douguet 
detonation.  Subjecting  the  differential  forms  of  the  (e-p-v)  equation 
of  state  and  the  Hugoniot  Eq.  (I-F-5)  to  the  condition  that  de  is  a 
perfect  differential  gives  the  following  equation  for  the  slope  of  a 
Hugoniot  curve  centered  at  (Po,v0): 
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FIGURE  1-1  TANGENCY  OF  THE  RAYLEIGH  LINE  TO  THE  PRODUCT 

HUGONIOT  CURVE  AT  THE  CHAPMAN-JOUGUET  <CJ}  POINT 


l-F-5 


Oe/dv)^  + J(P!  + to) 
{ &e/3p>  ^ - i<v0  - V,,) 


(I-F-7) 


Equation  (I-F-7)  gives  Oe/ap)  - |<v0  - vx)  / 0 as  a necessary  coodi- 

V 

tion  for  the  slope  of  the  Hugoniot  curve  to  be  finite.  Elimination  of 

(Se/Sv)  and  (pi  to)  from  Eq.  (I-F-7)  with  the  identity 
P 


a . „ P + (ae/&v) 

-'ey  _ _ /apN  p 

Vv/  Vav. / <ae/ap)v 


(I-F-8) 


and  the  equations  for  the  Rayleigh  line  written  as 

<Pi  “ Pq)  ~ t(D  - ui)/Vi}3  <va  - Vi) 
gives  the  equation 

(c1/v1)3<ae/ap)  - [(D  - Uii/vj]3  J(v0  ' vx) 

op  _ _ y ( 

dv  <ae/ap)y  - §(v0  - vx) 

Since  (de/dp)  - $<v0  - vx)  t 0,  combining  Eq.  (I-F-9)  with  the 
Chapoan-Jouguet  condition  Eq.  (I-F-6) , and  Eq.  (I-F-8)  gives  the 


( I -F  -9) 


equation 


*£  . 'Sx.)*  = (*E) 

dv  \v,/  \5v/ 

* c 


Thus,  the  Rayleigh  line  through  (to ,Vo)  that  satisfies  the  Chapman - 
Jouguet  condition  is  tangent  to  the  Hugoniot  curve  centered  at  (to*vo^ 
and  is  tangent  also  to  the  isentrope  passing  through  the  Chepraan->3ouguet 
point.  In  other  words,  the  Hugoniot  curve  has  a first-order  point  of 
contact  with  the  isentrope  passing  through  the  Chapmen -Jouguet  point. 

It  follows  from  the  tangency  condition  that  the  Chapman-Jouguet  detona- 
tion wave  propagates  with  a minimum  velocity.  The  Chapman-Jouguet 
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hypothesis  lends  to  the  conclusion  that  all  unsupported  reactive  shock 
discontinuities  satisfying  the  conservation  laws  of  mass,  momentum,  and 
energy  are  unstable  with  respect  to  the  discontinuity  propagating  at  the 
lowest  velocity, 

2.  States  Attainable  from  the  Chapraan-Jouguet  State 

Other  processes  resulting  from  the  interaction  of  a detonation 
wave  with  its  surroundings  are  important  for  understanding  the  behavior 
of  explosives.  A particular  interest  in  states  attained  by  the  reflection 
of  a shock  on  the  one  hand,  and  of  a rarefaction  wave  on  the  other  hand, 
into  a detonation  leads,  respectively,  to  the  construction  of  the  Hugoniot 
curve  of  detonation  products  centered  at  the  Chapman -Jouguet  point  and  to 
the  isentrope  of  detonation  products  through  the  Chapman -Jouguet  point. 
This  Kugoniot  curve  centered  at  the  Chapnan-douguet  point  specifies  the 
maximum  pressure  in  detonations  reflected  from  metals  with  increasing 
shock  Impedance;  the  isentrope  specifies  the  succession  of  states 
attained  in  detonation  products  expanding  into  the  atmosphere. 

3.  Equilibrium  Chapman-Jouguet  Calculations  with  an  Arbitrary 
Equation  of  State 

Our  treatment  of  CJ  states  will  be  based  on  more  general 
expressions  of  the  jump  conditions  than  those  given  in  Section  I-D-lb. 
Let  w denote  the  mass  velocity  with  respect  to  the  reactive  shock 
discontinuity,  and  let  the  subscripts  0 and  x again  denote  conditions 
at  the  bottom  and  the  top  of  this  discontinuity.  Then  the  Rankine- 
Hugoniot  jump  conditions  expressing  the  balance  of  mass,  momentum,  and 
energy  across  the  discontinuity  can  be  written  as: 


Plwl  = Pc*  w0 

pi  + 0i  wi S = Po  + 0bwo2 


hi 


+ £ W;3 


ho  + ! 


I-F-7 


(I-F-ll) 
(I-F-12) 
( I -F  -13) 
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The  Chmpman-Jouguet  condition  is  expressed  by  the  equation 


W!3  = Cl*  (I-F-14) 

with  the  sound  speed  c defined  by  the  equations 

ca  = -va(dp/dv)g  = xpv  (I-F-15) 

and  the  adiabatic  index  x related  to  the  thermodynamic  coefficients 
a and  0 by  the  following  expression  given  previously  in  Section  I-B-4e 


C 

v'  5p\  p V'  &p\  _ 0!  + 1 

pV  Sv/  C p\  Sv/  g 

s v T 


(I  -P-16) 


Combination  of  Eqs„  tl-F-ll)  and  (I-F-12)  leads  to  the  equation 


£2.  a JSL_  i'2ft  _ j \ 

Pi  Pi  vl'vl  } 


(I-F-17) 


and  combination  of  Eqs.  (I-F-ll)  through  (I-F-13)  leads  to  the  Hugoniot 
equation  (I-F-5)  , 1 - ec  = JCpj  + Po)(v0  - vj)  , relating  the  thermo- 

dynamic states  connected  by  the  shock  discontinuity.  Equations  (l-F-14), 
(I-F-15) , and  (I-F-17)  can  then  be  combined  to  give  the  following 
expression 


■ v»u7(l  + 


(I-F-18) 


for  the  CJ  condition  in  terms  of  thermodynamic  quantities  only.  Chapman- 
Jouguet  states  can  therefore  be  considered  as  thermodynamic  states 
defined  by  the  Eqs.  (I-F-5)  and  (I-F-18).  Since  TIGER  was  designed  to 
perform  thermodynamic  calculations,  these  equations  are  used  to  compute 
CJ  points  in  the  cod?,  but  an  iterative  procedure  must  be  used  to  solve 
them  because  ej  and  cannot  be  evaluated  at  a (pj point 

unless  the  composition  is  known. 


The  equations  used  to  step  forward  from  one  equilibrium  thermody- 
namic point  calculation  to  the  next  in  the  iterative  scheme  will  now  b.- 
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derived.  Let  superscript  prime  denote  a point  where  the  thermodynamic 
state  has  been  calculated  and  the  values  of  £e1/,  Pj ' , v, r , Tj7, 

Ctl  ' , fa  ' } are  known.  Then  Eq.  (I-F-18)  and  an  approximation  to 

the  Hugoniot  equation  are  used  with  these  values  to  generate  the  two 
state  variables  that  are  required  to  perform  the  next  thermodynamic 
point  calculation.  Let  Ap  = Pi  - p^ / and  Av  = then  combi- 

nation of  Eq.  I-F-5)  with  the  Taylor  series  expansion  for  et  about 
the  point  (p*  ' ,Vi  *)  , 

- ei'  + Oi  'f>i  'iv  + 'ip  (I  -F-1S 


gi  ves  the  equation 

atej'  - ec)  - (Pi*  + Po)(v0  - vx  ')  - [v0  - v1  ‘ - 2/3i'v1,'JAp 

(I  -F -2C 

“ [ Po  + Pi  * + 20^  #Pi]av  - apAv 

which  on  rearrangement  gives  the  approximation  to  the  Hugoniot  equation 
that  is  used  in  the  code 


Ap  = 


Cei  / - ec)  - (Pi'  + Po)  ( Vp  - v^)  f 2 
v0  - Vi  - 2vx  ' 


Equation  (I-F-21)  is  used  to  step  forward  in  pressure,  and  Eq.  (I-F-18) 
written  as 


"1  BlN  ,1 

Vl  v=iJ ' 70 ' J 


is  used  to  step  forward  in  volume  in  performing  the  sequence  of 
thermodynamic  calculations  to  determine  the  CJ  point , 

The  procedure  for  calculating  the  CJ  point  on  the  detonation  branch 
of  an  equilibrium  products  Hugoniot  curve  of  a condensed  explosive  is 
described  below.  The  case  when  no  problems  are  encountered  in  the 
calculations  is  discussed  first,  and  then  the  types  of  problems  that 
arise  in  such  calculations. 
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derived.  Let  superscript  prime  denote  a point  where  the  thermodynamic 
state  has  been  calculated  and  the  values  of  Pi  * * vi ' > Ti ' , 

Cti  ‘ < fa  ' . •••}  ate  known.  Then  Eq . (I-F-18)  and  an  approximation  to 
the  Hugoniot  equation  are  used  with  these  values  to  generate  the  two 
state  variables  that  are  required  to  perform  the  next  thermodynamic 
point  calculation.  Let  ip  = pA  pj  ' and  Av  = vx  - ' , then  combi- 

nation of  Eq.  I-F-5)  with  the  Taylor  series  expansion  for  e!  about 
the  point  (px  '.Vi  ')  , 

el  - ei  * + ®i  *Pi  V + fa  ,v1  'ap  (I-F-19) 

gives  the  equation 

2(6!'  - e0)  - (Pi  ' + p0)(v0  - vx  ')  = Cv0  - vr  / - 2^ 'vx 'lup 

( I -F  -20) 

~ [Po  + Pi'  + 2«i/Pi]iv  - ipiv 


which  on  rearrangement  gives  the  approximation  to  the  Hugoniot  equation 
that  is  used  in  the  code 


6p  „ ?ie.i.L:  ,e.a>  r.J.,Pi_....t-.ps>-<.yq.  - yx>.+-2pA_AY. 

V0  - vx  - 2V!  01  ' 


(I  -F  -21) 


Equation  (I-F-21)  is  used  to  step  forward  in  pressure,  and  Eq.  (I-F-18) 
written  as 


vi  * (**§{)*  0 


( I -F  -22) 


is  used  to  step  forward  in  volume  in  performing  the  sequence  of 
thermodynamic  calculations  to  determine  the  CJ  point. 

The  procedure  for  calculating  the  CJ  point  on  the  detonation  branch 
of  an  equilibrium  products  Hugoniot  curve  of  a condensed  explosive  is 
described  below.  The  case  when  no  problems  are  encountered  in  the 
calculations  is  discussed  first,  and  then  the  types  of  problems  that 
arise  in  such  calculations. 
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Determination  of  the  CJ  point  involves  the  sequential  performance  of 

an  (e,v)  point  calculation,  a (T,v)  point  calculation,  and  a series  of 

<p,v)  point  calculations.  The  (e,v)  calculation  is  performed  first 

with  ej  = e0  and  vt  = v0  to  obtain  the  constant  volume  explosion 

point  (e0>  v0  , p / , T ' , « / , A ' , . . .)  as  a first  approximation  to  the 
e e e re 

CJ  point.  The  (T,v)  calculation  is  then  performed  to  obtain  a better 

approximation  to  it.  Originally  the  (T,v)  calculation  was  performed 

with  Ti  = T ' and  v<  = v0x  V(x  ' + 1) , obtained  by  ignoring  the 
e e e 

pressure  term  pp/p  ' in  Eq.  (I-F-22) , but  experience  has  shown  that  a 
e 

better  approximation  to  the  CJ  state  is  obtained  by  choosing  Ti  = T ' 

e 

and  va  2 0 . 85  v0  . 

The  scries  of  (p[V>  calculations  are  performed  next  to  determine 
the  CJ  point.  The  first  (p,v)  point  calculation  is  performed  with 
values  of  px  and  vj  generated  fromEq.  (I-F-22)  and  with  the  values 
of  the  thermodynamic  variables  obtained  in  the  (T,v)  calculation. 
Specifically,  Eq . (I-F-22)  is  first  used  to  calculate  vx , and  <hen 
Eq.  (I-F-21)  is  used  to  calculate  the  corresponding  value  of  pa . The 
state  variables  calculated  at  the  first  (p,v)  point  are  then  substituted 
into  Eqs.  (I-F-22)  and  (I-^?-21)  to  generate  the  values  of  px  and  Vj 
for  the  second  (p,v)  calculation,  and  the  process  is  contained  until 

•43  8 

the  CJ  point  is  obtained.  When  l^v/vjJ  < 10  and  |Ap/pl|  £ 10 
between  successive  calculations,  it  is  assumed  that  the  CJ  point  has  been 
calculated  successfully.  It  is  clear  that  the  method  for  calculating  the 
CJ  point  on  the  detonation  branch  of  the  Hugoniot  can  easily  be  modified 
to  calculate  the  CJ  point  on  the  deflagration  branch  of  the  Hugoniot 
curve . 

The  original  iterative  scheme  for  determining  CJ  points  by  perform- 
ing a series  of  (p,v)  point  calculations  was  developed  for  gases  by 
S.  R.  Brinkley,  Jr.  A series  of  computational  problems  encountered  with 
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CJ  calculations  when  the  code  was  being  developed,  however,  led  to  the 
conclusion  that  this  scheme  was  unsatisfactory  for  condensed  explosives. 

L.  B.  Seely  solved  these  problems  by  introducing  the  (T,v)  point 
calculation  to  obtain  a better  approximation  to  the  CJ  point  in  the 
early  stages  of  the  iterative  procedure. 

Calculation  of  the  CJ  point  becomes  a problem  when  the  code  cannot 
compute  the  constant  volume  explosion  state  or  the  state  at  one  of  the 
(p,v)  points  in  the  series  of  (p,v)  calculations.  The  constant  volume 
explosive  state  could  not  be  computed  for  nonideal  explosive  compositions 
because  it  lay  in  the  mixed  phase  region  of  one  of  the  condensed  con- 
stituents. No  difficulties  have  been  found  with  these  compositions, 
however,  since  the  revision  of  the  routine  for  handling  systems  contain- 
ing multiple  condensed  species  capable  of  changing  phase. 

Difficulties  in  the  (p,v)  points  calculations  were  found  to  arise 
for  the  following  reasons:  (1)  the  r's  do  not  converge,  (2)  the  tem- 

perature is  out  of  range,  (3)  the  n.'s  do  not  converge,  (4)  the  code 
cannot  find  a set  of  ondensed  components,  and  (5)  the  set  of  equations, 

F — 0 with  3=1,  . c + 2,  has  no  solution  because  of  a local 
minimum.  In  the  first  three  cases,  the  step  size  Av  taken  to  reach 
the  incalculable  point  is  halved  until  convergence  is  obtained,  and  then 
the  series  of  (p,v)  point  calculations  is  continued.  The  last  two 
cases  arise  when  the  values  p and  v chosen  for  the  (p,v)  calcula- 
tion are  incompatible  with  the  equilibrium  states  of  the  system.  In 
either  case,  a (p,T)  calculation  is  performed  to  obtain  an  equilibrium 
(p,v)  point,  and  the  iterative  procedure  with  (p,v)  points  is  continued. 
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4.  Construction  of  the  Hugoniot  Curve  for  an  Arbitrary  Equation  of  State 


The  iterative  procedure  for  constructing  a point  on  the  Hugoi  lot 
curve  by  performing  a series  of  (p,v)  calculations  when  either  the 
volume  or  pressure  is  specified  is  based  on  Eq,  (HP-20)  . 


Consider  first  the  case  when  the  thermodynamic  state  at  a (p^v') 
point  is  known,  and  the  Hugoniot  point  at  a specified  volume  vx  is 
required.  Equation  (I-F-21)  is  used  to  generate  the  first  approximation 
to  the  Hugoniot  point.  Then  the  equation 


Ap 


2(e!  ' - e0)  - (p!  ' + pp)  (v0  - 
v0  - Ml  + ' 


(I-F-23) 


obtained  by  setting  Av  = 0 in  Eq.  (I-F-21)  is  used  to  generate  sub- 
sequent approximations.  The  iteration  process  is  started  by  substituting 
the  known  values  of  Av,  e‘,  p ‘ , o! , and  fi'  into  Eq.  (I-F-21)  to  obtain 
the  value  of  pj  for  the  next  (Pi ,vx)  calculation.  This  procedure  is 

continued,  using  Eq.  (I-F-23) , until  the  Hugoniot  point  is  obtained.  When 
-e 

lAP/Pii  10  between  successive  calculations,  it  is  assumed  that  the 
Hugoniot  point  has  been  calculated  successfully.  The  Hugoniot  curve 
passing  through  the  point  calculated  at  Vj  can  be  readily  constructed 
by  performing  a series  of  such  calculations  for  different  values  of  v. 

But  then  it  is  convenient  to  use  the  equation 


Ap 


t Pi  *(1  + 20k  *)  + Pol  Av 
v0  - Vl(l  + 2 & ') 


(I  -F  -24) 


rather  than  Eq.  (I-F-22)  to  find  the  value  of  pt  for  the  first  (p,v) 
calculation  because  the  original  state  that  is  known  lies  on  the  Hugoniot 
curve . 
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Now  consider  the  case  when  the  thermodynamic  state  at  a (p',v') 
point  is  known  and  the  Hugoniot  point  at  a specified  pressure  pi  is 
required,  A similar  procedure  is  followed  with  fixed  and  values 

of  evaluated,  using  the  equations  for  Av  obtained  from  (I-F-20). 

The  value  of  vx  for  the  first  (Pi calculation  is  obtained  from 
the  equation 


Av  = 


+ Pto)  (Vo  - Vi  *)  - 2(ei  ' - eQ  + ft  'vi  'ap) 
PO  + Pi  + 20^  p*  ' 


C I -F  -25) 


corresponding  to  Eq . (I-F-21),  and  values  of  v for  subsequent  (pi  ,v) 
calculations  are  obtained  from  the  equation 


Av  = 


+ - Vi  ')  - 2(e,  ' - e0) 

Po  + Pi  < 1 + 20!  ') 


(I  -F  -26) 


corresponding  to  Eq.  (I-F-23)  . When  the  known  point  lies  on  the  Hugoniot 
curve,  the  value  of  v for  the  first  (p,  ,v)  calculation  is  obtained 
from  the  equation 


Av  = 


Cvo  - Vi  'U  + 2ft  ')]/ 
Pc  ~+  PiU  + 20(i 


( I -F  -27) 


corresponding  to  Eq . (l-F-24).  As  in  the  previous  case,  Hugoniot  points 
are  assumed  to  be  determined  when  | Av/vj | ^ 10  between  successive 
(Pi , .x)  calculations. 

Construction  of  the  Hugoniot  curve  becomes  a problem  when  a (p,v) 
point  is  incalculable  because  of  discontinuities  associated  with  a con- 
densed species  disappearing  from  the  system,  in  one  case,  the  incalcul- 
ability  arises  because  the  code  cannot  find  a set  of  condensed  components; 
in  the  other,  because  there  is  no  solution  to  the  set  of  F^  = 0 equations 
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for  the  specified  values  of  p and  v.  In  either  case,  the  values  of 
p and  v chosen  for  the  (p,v)  calculation  are  incompatible  with  the 
equilibrium  states  of  the  system,  and  the  problem  is  treated  in  the  same 
way  as  in  the  CJ  routine.  A (p,T)  calculation  is  performed  to  obtain 
an  equilibrium  (p,v)  point,  and  then  the  iterative  procedure  with 
(p,v)  {joints  is  continued. 

5.  Construction  of  an  Isentrope  for  an  Arbitrary  Equation  of  State 

An  isentrope  is  readily  constructed  by  performing  either  a series 
of  (s,v)  or  a series  of  (s,p)  point  calculations  with  a constant 
value  of  s.  The  CJ  isentrope  is  obtained  by  starting  the  calculations 
at  the  CJ  point, 

6 . Frozen  Chapman-Jouguet  Calculations  with  an  Arbitrary  Equation  of 
State 

Frozen  CJ  states  are  readily  calculated  in  TIGER  by  performing  non- 
equilibrium point  calculations  in  the  CJ  routine  with  the  partial  freeze 
options  described  in  Section  I -C  of  the  TIGER  documentation.  The  frozen 
CJ  states  calculated  are  therefore  those  in  mechanical  and  thermal  equi- 
librium, with  part  of  the  composition  frozen  and  part  of  it  in  chemical 
equilibrium.  The  mole  numbers  of  the  frozen  species  are  fixed,  but 
frozen  condensed  species  are  allowed  to  change  phases  according  to  their 
equations  of  state  as  in  equilibrium  calculations.  Frozen  calculations 
can  easily  be  performed  without  phase  changes  by  removing  the  equation 
of  state  of  one  of  the  phases  from  the  library. 

The  SPECIAL  routine  for  computing  nonequilibrium  CJ  states  with  two 
frozen  constituents  was  developed  to  provide  Picatinny  Arsenal  with  a 
capability  to  investigate  nonideal  explosives.  Ix  was  developed  primarily 
for  computing  frozen  CJ  states  with  part  of  the  explosive  composition 
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unreacted,  but  tt  can  obviously  be  used  for  computing  CJ  states  when  the 

explosive  is  completely  decomposed  and  two  of  the  reaction  products  are 

frozen.  The  reaction  coordinates  Ax  end  X.2  were  introduced  to  describe 

the  partially  frozen  nonequt librium  states.  The  coordinates  A,(i  = 1,2) 

were  chosen  to  satisfy  the  conditions  0 ■£.  A <1  by  setting  A = a /m  , 

i i l x 

where  a.  denotes  the  mass  of  species  i frozen  in  1 kgm  of  mixture, 

and  m denotes  the  largest  value  of  a that  will  be  used  in  calcu- 
i i 

lations  with  the  particular  composition.  Since  in.  can  be  chosen  to  be 
larger  or  smaller  than  the  equilibrium  mass  of  species  l , frozen  CJ 
states  containing  more  or  less  of  species  i than  the  equilibrium  CJ 
state  can  be  considered. 

SPl'CIAL  was  written  to  provide  a complete  thermodynamic  description 
of  frozen  CJ  states  with  different  values  of  Ax  and  Xs . Included  in 
this  description  are  the  parameters  listed  in  the  equilibrium  CJ  routine, 
the  rate  of  entropy  production,  two  additional  quantities  (Ox  and  03 
(which  can  be  preselected) , the  first  and  second  derivatives  of  these 
quantities  with  respect  to  volume  along  a frozen  Hugoniot  curve,  and 
the  first  derivatives  along  the  loci  of  frozen  CJ  states.  These  first 
and  second  directional  derivatives  are  approximated  by  differencing 
techniques  when  their  explicit  expressions  are  not  available. 
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